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SUMMARY 


A computational method for simulating the aerodynamics of wing-fuselage config- 
urations at transonic speeds has been developed. The finite difference scheme is char- 
acterized by a multiple embedded mesh system coupled with a modified or extended 
small -disturbance flow equation. This approach permits a high degree of computational 
resolution in addition to coordinate system flexibility for treating complex realistic air- 
craft shapes. To augment the analysis method and permit applications to a wide range 
of practical engineering design problems, an arbitrary fuselage geometiy modeling 
system has been incorporated as well as methodology for computing wing viscous effects 
Configuration drag is broken down into its friction, wave, and lift -induced components. 
Typical computed results for isolated bodies, isolated wings, and wing-body combina- 
tions are presented. The results are correlated with experimental data. A computer 
code which employs this methodology is described in the Appendix of this report. 




INTRODUCTION 


The inherent complexity of transonic flows renders a purely experimental 
approach to solving aircraft design problems impractical in many cases. Mixed 
subsonic and supersonic regions, shock waves, and complicated viscous effects 
characterize this flow field. Numerical methods for simulating transonic flows are 
being developed to reduce the design cost and improve the performance of high speed 
aircraft. 

Murman and Cole ^ developed the first practical technique for computing two- 

dimensional inviscid transonic flows. This scheme utilized mixed (upwind and central) 

finite differences which were applied to the small perturbation equation in Cartesian 

coordinates. The equation was solved iteratively by numerical relaxation. The success 

(2 - 4 > 

of this technique made it possible for other investigators ' ’ to develop methodology 

incorporating the full potential equation and exact boundary conditions. Typically, 
coordinate system transformations now’ provide a grid system which is closely aligned 
with the airfoil surface. The mapped coordinates simplify application of the exact 
boundary conditions which are consistent with the full potential equation method. One 
advantage of this approach is that details of the high -gradient region near the airfoil 
leading edge are accurately predicted. 

It is interesting to note that three-dimensional numerical methods for computing 

transonic flows about wings and wing-body combinations have evolved in the same 

( 5 ) 

manner that airfoil methods evolved. The early work of Bailey and Steger with 
the transonic small disturbance equation provided the foundation for a number of three- 
dimensional schemes now in existence ^ ^ . Full potential equation codes have also 

been developed for treating wings and simple wing -body shapes ^ Unfortunately, 

the geometry-fitting coordinate system which is very advantageous when applied to 
simple airfoil shapes becomes increasingly difficult to implement as geometries become 
more complex. Constraints on computational resolution and geometry/coordinate 
system flexibility arise in three-dimensional applications. 

A high degree of numerical resolution will be necessary if flow details about 
complex shapes are required. Providing this resolution for both a wing and fuselage 
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simultaneously presents a difficult problem. In addition, computational mesh flexi- 
bility is necessary to permit the treatment of truly arbit rary shapes. This must all 
be accomplished with sufficient computing efficiency to insure a cost-effective 
approach to solving engineering design problems. These requirements place a 
severe strain on the conventional single continuous grid/coordinate transformation 
approach which characterizes existing numerical methods. 

This grid embedding approach provides a means for modeling and analyzing 
configuration components (wings and bodies) within local detailed fine grid regions. 
Beyond these fine grid regions, a global crude mesh treats the flow field. Fine mesh 
computations are performed only in a region very close to the configuration where 
gradients are large and details are important. This permits a high degree of computa- 
tional resolution and, at the same time, reduces computer resource requirements. 
Perhaps more important, mesh flexibility for treating complex shapes is available 
since global transformations are not required. 

The present effort provides methodology for treating the more complex realistic- 
aircraft shapes that occur in pratical applications. The method uses a modified or 
extended small disturbance flow equation. With full potential equation solutions in 
hand, this would at first appear to be a step in the wrong direction. However, the 
simplicity of the planar boundary conditions which accompany the small disturbance 
equation, permits the extensive use of Cartesian coordinates; this, in turn, 
facilitates the implementation of a grid embedding scheme. 
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NOMENCLATURE 


SYMBOLS 

C 

C AV 

C D 

C D F 

C Di 

Cd body 
C DP (b) 

C Dp (w) 
Cn 

d wing 

Cn 

U WAVE 

° f BODY 

Cf wing 

C L 

C LBODY 

Cl wing 

C M 

Cm body 

C MWING 

C R 

c “ b 


Infinite sheared wing chord 
Average wing chord 
Total configuration drag coefficient 
Configuration friction drag coefficient 
Configuration lift-induced drag coefficient 
Body drag coefficient 

Body pressure drag (lift-induced and wave drag) 

Wing pressure drag (lift-induced and wave drag) 

Wing drag coefficient 

Configuration wave drag coefficient 

Body skin-friction drag coefficient 

Wing friction drag coefficient 

Total configuration lift coefficient 

Body lift coefficient 

Wing lift coefficient 

Total configuration moment coefficient 

Body moment coefficient 

Wing moment coefficient 

Wing root chord at configuration centerline (<£) 

Wing tip chord 

Body cross-section pressure drag coefficient (lift-induced 
+ wave drag) 



SYMBOLS 


NOMENCLATURE (Continued) 


C d 

U) 



W 


Cm, 

b 

Cm 

u) 


-'f n 


w(loc) 


Cp 


Cp* 

Cp C 2> C 3 , B lf B 2 
F 


L 


M. A.C. 


M 


00 


Wing section pressure drag coefficient (lift-induced 
+ wave drag) 

Wing section skin-friction drag coefficient 

Body cross-section lift coefficient 

Wing section lift coefficient 

Body cross-section moment coefficient 

Wing section moment coefficient (taken about reference 
position) 

Wing section moment coefficient (taken about local 
quarter-chord position) 

Pressure coefficient 

Critical (sonic) pressure coefficient 

Crude grid scaling constants 

Body geometry function 

Body length 

Mean aerodynamic chord 
Freestream Mach number 


N x ,N y , N z 


R 


Re 

S 


S BASE 

S PB 

SBwET 


Surface normal direction cosines 
Body surface distance 
Reynolds number 
Arc length 

Body maximum cross-sectional area 
Body projected area (x-y plane) 

Body wetted area 


S 


cs 


Cross-sectional area 
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NOMENCLATURE (Continued) 


SYMBOLS 

S EXP 

S W 

X C/4 

X REF 

b 


E 

f 


u,v,w 


x,y,z 

Ax, Ay, Az 
00 

r 

A c/2 

a 

7 

^ * 


X 

£ , , r 

p 

T 


Exposed wing area (for wing-body combinations) 

Wing area 

Wing section quarter-chord position 
Position about which moments are computed 

Wing span 

Oswald efficiency factor (lift induced drag) 

Wing geometry function 
Perturbation velocities 
Physical coordinates 
Mesh spacing 
Infinity 
Circulation 

Wing mid-chord sweep angle 
Angle -of-attack 
Specific heat ratio 

Wing boundary layer displacement thickness 
Wing span position ( ij = 2y/b) 

Momentum thickness 
Wing taper ratio C T /C R 

Computational coordinates (symbolic for all grid systems) 

Density 
Shear stress 

Perturbation velocity potential 


SYMBOLS 


NOMENCLATURE (Concluded) 


U) 

F 

SUBSCRIPTS 

GLE 

GTE 

TE 

W 

b 

c 

e 

1 

loc 

u 



SUPERSCRIPTS 


Subsonic relaxation factor 


Body lineness ratio 


Grid leading edge 
Grid trailing edge 

\\ ing section trailing edge 
Wing 

Body surface 
Computational surface 
Edge condition 


Local wing section fine grid 


Lower surface 


Local 


Upper surface 


Mesh indices 


Partial derivatives 


Dummy value 

Infinite sheared wing coordinate 
Wing/Wake upper surface 
Wing/Wake lower surface 
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COMPUTATIONAL METHOD 


Flow Equation 

Early three-dimensional computations using the classical transonic small- 
disturbance equation indicated that shock waves with appreciable sweep (>15 5 ) relative 
to the freestream direction could not be resolved. As a result, Lomax, Bailey and 
Ballhaus (17) proposed the retention of additional terms in the flow governing equation. 
Their derivation was based on simple sweep theory concepts. 

The present method employs a flow equation constructed by adding three addi- 
tional terms to the classical transonic small -disturbance equation. The flow equation is 


[1 ~ M^,— (7 + 1) M^x - M °°^x^xx 2M ^y^xy 
+ [1 - (7 - 1 ) M^ x ] * yy + ^zz = 0 


All of these terms can be found in the full potential equation (see ref 18, page 204). 

The * * and * * terms , typically referred to as crossflow terms , permit the 
resolution of sho^k waves with large sweep angles. A graphic example of the effects 
of these terms on computed pressure distributions and shock wave pattern can be seen 
in Figure 9 of ref 16. The *** xx term has been retained to provide a better approxi- 
mation of the critical velocity at which the full potential equation changes type (elliptic 
to hyperbolic). Difference approximations are applied to equation (1) as it is written. 
Empirical corrections (19 > and similarity variables are not employed. The method 
could be classified as modified or extended small disturbance in character. 

Pressure coefficients on wing surfaces are computed using the following equation. 

Cd = — f2y>„ + (1 M 00 )y> x + 'Pyl (2) 


To simplify velocity computations on the non- 
is used. Note that the* and* components as 
eluded in the present version of the method. 


planar body surface, a simplified equation 
listed in ref 18, page 206 are not in- 


Cp = — [2<p x + (1 M2) * x ] 


(3) 


9 



Computational Space and Grid Systems 

The computational space used in the present method is illustrated in Figure 1. 
This space is filled with a relatively crude Cartesian mesh. Instead of adopting a far- 
field solutlon for the grid outer boundaries, the original x, y, z region is stretched to 
a ^ T region in which the boundaries correspond to infinity. The flow field poten- 
tial is set to zero on all bounding planes except the dowmstream plane for which the 
following equation is solved. 


^yy + ^zz ® 


The following conditions are enforced at the symmetry plane . 


^ y = ° 
^xy “ 0 


(4) 


(5a) 


(5b) 

Carlson (20) has noted that caution must be exercised in selecting a grid stretch- 
ing function. For this reason, simple tangent functions are employed. Space in the 

x -direction is broken up into three regions (see Figure 2). The grid stretching takes 
the following form: 

Regions I and III X - Xj ♦ B, TAN <§ (f - {l)( ♦ B 2 TAN ( | (f ^ )3| ^ 

Region II X = C.f 

(6b) 

Note that an evenly spaced grid spans the wing planform tor wing and wing-body cases* 
For isolated bodies , the entire body is covered by Region II. The computational space 
is divided into 50 cells or segments in the x-direction. Of these. 38 cells are con- 
stramed to fall within Region II. The constants Cj and B arc computed by requiring 


X » Xj at $ - 


(6 c 


Tim author would like to thank Ed Waggoner (Vought Corporation, for demonstrating 

improved performance of this x-grid system over that originally employed 
(see ref 16). J 
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and 


atX = X 1 

Wl \ d ^ /II (6d) 

The constant B 9 is determined by requiring that the next -to -last physical x-grid point 
be stretched several configuration lengths upstream and downstream. As a xesult, the 
finite computational region: 


-4. 167 <$ <4.167 


(6e) 


corresponds to the infinite physical region: 

— o° < X < 00 (6f ) 

Grid stretching in the y-direction takes the following form: 


TANH ' 1 (77) 

c 2 (7a) 


The hyperbolic tangent saves considerable computing resources >vhen compared 
to the milder tangent function. The use of this exponential function in the spanwise 
direction has resulted in no apparent problems. The constant is adjusted to 
provide 18 crude grid span lines along the wing out of the total of 26. In addition, it 
is required that the wing tip lie exactly at the midpoint lietween two y— grid lines. This 
central tip location provides the best approximation of the wing tip aerodynamics. The 
finite computational region: 

0 < V < 1 (7b) 


now- corresponds to the infinite physical region: 

0 < y < 00 (7 C ) 


The grid stretching function in the z-direetion is: 

Z = c 3 TAN (Jf ) (8a) 

Thirty-one grid lines are employed with the constant C chosen to provide sufficient 
mesh stretching so as not to inhibit or suppress supersonic flow 1 development. As a 
result, the finite computational region: 

-1 < ? < 1 (8b) 
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now corresponds to the infinite physical region: 


— 00 <C Z <C 00 


(8c) 


The stretching adds transform derivatives to the crude grid governing equation. 
The following substitutions are made in equation (1). 


£ x 

(9a) 

^xx “ + £xx 

(9b) 

II 

<< 

(9c) 

^yy ~ ^vv ^y^ 2 + ^yy 

(9d) 

^xy ~ ^£77 ^x^y 

(9e) 

^zz = % (f z ) 2 + ^$zz 

(9f) 


This results in a total of 41, 106 grid points representing the global crude grid. 
The primary purpose of this mesh system is the proper representation of both configu- 
ration effects in the far-field and far-field conditions on the configuration. However, 
the crude mesh also serves as a communication link for the embedded fine grid arrays 
to be described. 


Individual fine grid arrays are constructed for the wing and body. These second- 
ary mesh systems serve two purposes. First, detailed computations are performed 
only in a region very close to the wing or body where gradients are large and details 
are important. The resulting numerical efficiency permits a very dense computational 
mesh, a tenefit in both the resolution of shock waves and the calculation of configu- 
ration forces and moments. Second, the embedded mesh systems are independent 
and optimized for a particular geometric component (wing or body shape). The system 
is not constrained by a single geometry-fitting transformation. This will facilitate 
future applications to configurations with multiple wing (wing, winglets , pylons, 
canards) surfaces and multiple tody (fuselage, pods, engines) surfaces. 


13 


Wing fine grid arrays are set up at each position where a crude y-grid line 
cuts the wing surface. This results in 18 fine arrays for isolated wings and somewhat 
less (depending on fuselage span extent) for wing -body combinations. Figure 3 is a 
schematic which illustrates the wing section/grid system arrangement at a particular 
span station along the wing. 



Figure 3. Wing Section/Crude Grid/Fine Grid Boundary Arrangement in 
Physical X-Z Plane 
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The wing fine grid is evenly spaced in both the x and z-directions (see Figure 4). 
As a result, a shock wave will be resolved to an equivalent degree independent of its 
position along the wing chord. The wing section leading and trailing edge are positioned 
at the mid-point between two streamwise mesh points. The present formulation 


utilizes 100 mesh points along the wing chord. Fine grid boundaries are positioned at 

20 ' l0Cal Ch ° rd dlstance in front of each wing section leading edge and 10% behind each 
trailing edge. Similarly, boundaries are set at 30% of the wing average chord distance 
above each section and 10% below. These limits have been pre-set in the solution coding 
and were used for all computed results presented herein. They may, of course, be 

modified for special applications . The total number of field points in the wing grid 
system is approximately 60,750. 



Figure 4 Details of Fine Embedded Wing Grid 
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The wing fine grid system has been sketched in Figure j. Note that the embedded 
grid system is tapered as is the wing planform. A shearing transformation is intro- 
duced so that neighboring fine grid arrays can be used directly for cross -flow 
differencing. 


x = £ ( X GTE _ X GLE> + X GLE 
= TANH 1 ^) 
y ~ C 2 
z = ? 


(10a) 

(10b) 

(10c) 


As was the case for the crude grid stretching transformation, this shearing trans- 
formation introduces additional terms into the flow equation. The following substi- 
tutions are made in equation (1). 


^x = 

(11a) 

*xx = ( ^x> 2 

(lib) 

‘fiy = + *7] Vy 

(11c) 

'Pyy = £yy + hyy + <*y> 2 + 2 tyy + (V Y )2 

(lid) 

fi’xy = ^xy + ^x^y + Vx 

(lie) 

'fizz = % 

(Ilf) 


Only the term ^ which results from the ^ term (Equation lib) is upwind 

differenced at supersonic points (see finite difference approximations). AH other 

. ( 11 , 21 ) 

terms are centrally diiferenced 

The body fine grid system is constructed to completely encapsulate the body or 
fuselage shape. Figure 6 illustrates the body fine grid arrangement. Fine mesh bound 
aries are positioned at J0 r , body length in front of the nose and behind the tail, lor 
special applications, the fine mesh may lie concentrated on one particular region for a 
detailed analysis. Boundaries above , below and on the side are positioned at approx- 
imately one average body radius away from the body surface. 
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As was the case for the individual embedded wing grid systems, the body mesh 
is evenly spaced in all directions. Transformations are not used. As a result, the 
flow governing equation (1) is differenced directly. Details of this grid system can be 
seen in Figure 7. A total of 54,000 points are available in this grid system, which 
brings the total points in all grid systems combined to approximately 156,000. 

Overlap regions in which the flow is computed twice provide a means for inter- 
acting the different grid systems. This approach, which also employs combined 
Dirichlet and Neumann boundary conditions and dummy boundary surfaces, is described 
in the solution process section. 


BODY 

BOUNDARY 

POINTS 
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Finite Difference Approximations 


Finite difference approximations* are substituted for the derivatives in equation 

(1). At each mesh point, the local velocity is computed to determine whether the flow 

is subsonic (elliptic equation) or supersonic (hyperbolic equation). Equation (1) can be 
written 

T *xx + U *xy + V *Vy + *> zz = 0 

( 12 ) 

where 


T=1-m£ -(7 + 1) m2 , ^-111 M 2^2 

U = -2M 2 

V=1-( 7 -1)m2 ^ 


(13a) 

(13b) 

(13c) 


The elliptic equation results when 


DES = U 2 — 4 • T • V < 0 


(14) 


and the hyperbolic equation results for DES > 0. With the flow type determined, the 
appropriate central or upwind difference operator is selected. For subsonic points 

= co^bj.k 2 ^ j, k + ^i + 1, j, k 


<P- 


AX" 


(15) 


and for supersonic points 




XX 




AX" 


(16) 


Central difference operators are used for all other derivatives regardless of the 
local flow velocity. 


*x 


^i + l,j,k ^i — l,j. k 


2AX 


(17a) 


= ^i, j + 1, k 
^ 2AY 


'fh J ~ 1. k 


(17b) 


* Finite difference operators are listed for equally spaced grid system only. 
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(17c) 


^i.i.k + 1 *i,j,k-l 
2AZ 


gj- 1, i - 1, k ~ W~ 1, i + 1, k + 1,1 -- 1, k + + l,j+ 1, k 

^xy 4AXAY 

_ ^i, j — 1 . k~^i, j,k ^i, j, k + ^i, j + 1 , k 


<P- 


yy 


AY 2 


(I7d) 


(17e) 


*^zz 


i, J> 


1 — 2 ^i, j, k + j, k + 1 


AZ" 


(17f) 


The (+) superscript indicates that a new or updated value oi the flow field potential is 

used while no superscript indicates that the value is from a previous sweep or cycle. 

This procedure enhances the diagonal dominance of the numerical scheme improving 

(21 22 ) 

the stability of the solution ’ 

It should be noted that non-conservative difference operators have been selected. 
As a result, mass is not conserved across shock waves and the discrepancy or loss 
increases with increasing shock strength. However, shock waves are weakened in the 
physical flow field by interaction with the viscous boundary layer. It has been known 
for some time that nonconservative differencing fortuitously approximates the real flow- 
shock wave strength. Since the present viscous correction method does not have suf- 
ficient detail (approximately 41 ehordwise boundary layer points) to resolve the details 
of the true shock-boundary layer interaction, non-conservative difference operators are 
employed. This should result in a better simulation of the physical flow field. 

Wing and Body Boundary Conditions 

Surface boundary conditions are enforced on ali outer surfaces of the computa- 
tional space, wing and body surfaces and grid interface regions. Computational space 
outer boundary requirements were illustrated in Figure 1. ISoundary conditions at 
grid interface regions will be discussed in the solution process section. This 
section will primarily concentrate on the boundary conditions applied to the wing and 
body surfaces and the wake potential jump condition. 
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Boundary conditions are imposed by setting the value of a potential or its first 
derivative at a field point which represents the configuration surface. The wing and 
wake are approximated by a planar surface which extends from the wing leading 
edge to downstream infinity, while the body or fuselage is represented by a fixed cross 
sectional surface extending from upstream to downstream infinity. Corrections are 
applied to the body boundary conditions for the simulation of finite length bodies . 
Modeling is sufficiently flexible to permit the treatment of wings at varying height 
relative to the body (high-low-mid -wing). 

Special difference approximations are required at boundary points. The planar 
wing and wake are represented numerically by a grid surface of double valued mesh 
points (see Figure 8). For a wing surface defined by 

z = f (x,y) (18) 

the wing flow tangency condition is approximated by 


'Pz (x, y, o) - f x — a + 6* (ig 

where the slope of the boundary layer displacement thickness is added only for the 
inviscid/viscous interaction mode of operation. 

The wing upper and lower surface boundary conditions enter the solution 
formulation by way of the ^ term in equation (1). At the wing surface, this term 
can be written 


'Pzz 


_ 1 /^K + 1-^K\ 

AZ ) y AZ 



Note that the I and J subscripts have been dropped for convenience, 
the following relation into equation (20) 


( 20 ) 

By incorporating 


x X 2 AZ 


The wing boundary condition on the upper surface becomes 


( 21 ) 


<P- 


zz 


, +, 2 
(x,y,o ) = 7Tn 


AZ 


~/*K + 1 -*K\ 

l C C * \ 1 

LI az ) 

■ ( f *u-“ +s *)J 


( 22 ) 
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Similarly, the wing lower surface boundary condition becomes 


^zz ( x - y. o' 


2 1 

’ / - 1 \ 

-If -a- 6* \ 

AZ 

V AZ I 

\ x e x / 


( 23 ) 


At the end of each sweep of the flow field, the Kutta condition is enforced by 
calculating the circulation at the trailing edge of the wing section 


^ ^(X/jE, y» 0 + ) - ^(X-pg, y, 0 ) 


( 24 ) 
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The magnitude of the circulation is determined as part of the solution. The circula- 
tion jump condition is enforced along the wake between the wing section trailing edge 
and downstream infinity. The finite difference approximation for the wake upper sur- 
face is 


^zz 
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and for the wake lower surface 
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The use of the shearing transformation (Equation 10-a) for the wing embedded 
fine grid system complicates the imposition of symmetry conditions for isolated wings 
and wing root juncture conditions for wing -body combinations. The simple Cartesian 
(crude) grid symmetry condition (Equation 5a) becomes 


*y = n *y + = 0 

for isolated wings at the symmetry plane and 


(27) 


*y"*t* y + ^y =F x 


(28) 


for wing-body combinations in the wing root juncture region. Here, F x represents 
the slope of the fuselage or body at the wing root. 


Computations indicate that numerical instabilities will result if special attention 
is not given to the selection of difference approximations in this region. These 
difficulties result from the nature of the shearing transformation. To solve this 
problem, a plane of dummy mesh points is positioned across the symmetry plane 
(within the body surface for wing-body configurations) opposite the wing. These flow- 
field potentials are artificial in the sense that there is no physical flow field associated 
with them. They simply provide a side boundary of potentials which when used for 
differencing, produce the proper side condition given by Equation (27) or (28). 
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A special first order accurate one-sided difference operator is used to generate 

the dummy interior point potential values. Figure 9 is a schematic which illustrates 

the mesh point arrangement. For if grid lines which are swept back in the physical 

plane (£ <0), the following equation is used 

y 


D ,f y q y D 
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(29) 


and for v grid lines which are swept forward (f > 0) 
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(30) 


Note that the operator changes depending on whether the grid lines are swept forward 
or backward. In each case, the coefficient of the dummy potential at the point (I, J-l) 
is larger than the coefficients of other potentials in the difference equation (see ref 
23). This enhances the effective diagonal dominance of the system even though the 
dummy points are not directly relaxed in the conventional sense*. 



Figure 9 Wing Root Boundary Condition for Dummy Interior Points 


*Problem solution provided by J. South, NASA LRC. 
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Differencing at the wing tip is complicated by the fact that the fine mesh system 
does not extend beyond the wing tip. Unlike the conventional global transformation 
approach, the coordinate lines do not have to toe unswept or unsheared far from the 
wing. Provisions must be made, however, for properly ending the fine grid computation 
at the wing tip. For this reason, another temporary fine mesh is positioned just 
beyond the wing tip. Like its neighboring fine grid on the wing, this grid array is 
located at a crude mesh Y-line. Its extent in the streamw'ise and vertical directions 
is consistent with the fine grid system on the W'ing. Both the dummy plane beyond the 
wing root and the temporary fine grid beyond the wing tip are computed for each sweep 
of the array of fine grid structures. While the root dummy plane is computed using 
difference formulas, equations (29) and (30), the tip plane is simply filled using linear 
interpolation and potentials from the crude Cartesian grid. 

A computational surface represents the true body surface in the Cartesian grid 
systems. For axisymmetric bodies, the computational shape is an infinite cylinder 
with a radius equal to the true body maximum radius. For wing -body configurations 
the shape takes the form of the body cross-section in the wing-body juncture region. 
Figure 10 is a schematic illustrating the true and computational body surfaces. This 
approach proves to be more economical and stable than "fitting" the body boundary 
points to the true tody shape in the axial direction. This approximation is probably 
consistent with the planar witig/w'ake surface and the small disturbance approach in 
general. 

Body boundary conditions are imposed by computing appropriate flow field 
potentials at points near the tody surface (tody computational surface). The flow field 
is relaxed, holding these values fixed during each iteration cycle. No attempt is made 
to position mesh lines so points fall near the tody surface. The mesh density is 
sufficiently detailed to provide a good hands-off representation of the cross-sectional 
shape (see Figure 7). For a tody shape defined by 

F (x, y, z) = 0 ( 31 ) 

The small-disturbance boundary condition is 

F x + F y^y + F z^z = 0 n9 . 
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Following Bailey and Ballhaus three-point extrapolated differences are substituted 

for the flow velocities in Equation (32). The body/mesh arrangement for a typical body 
point is illustrated in Figure 11. 

Between the nose and tail of the body, there are sections where the true body 
surface and the computational body surface are not aligned. This is illustrated in 
Figure 10 and can be seen on a sectional basis in Figure 11. To ensure a proper rep- 
resentation of the true body shape, a correction must be applied to the boundary 
condition on the computational surface. The boundary condition correction will account 
for the displacement effect. Similarly, flow velocities computed on the computational 
surface must be corrected to yield velocities and pressures on the true body surface. 
Slender body theory provides the required modifications*. 


* Boundary condition displacement corrections provided by Professor Jack Werner, 

(2 4) 

Polytechnic Institute of New York ' . 
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Figure 11. Details of Body Surface Approximation 


The effects of body thickness and incidence lead to separate corrections. Body 
thickness effects can be approximately represented by a source distribution, the 
strength of which is proportional to the rate of change of body area. The boundary 
condition on the computational surface can be obtained by requiring that the net source 
strength across both surfaces be identical at each cross-sectional cut. 





( 33 ) 


As a result, the correct global effect of the true body shape on the remainder of the 
flow field will be obtained. This leads to 



( 34 ) 
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In the present method, cross -sectional shapes are taken to be circular or near 
circular in character. Equation (34) is approximated by 

“R 

c C b < 35 > 

A similar correction can be obtained for effects caused by flow incidence or 
angle -of -attack. In this case equivalent doublet strengths are enforced where doublet 
strength is proportional to the cross-sectional area S. 


«eff S c < X > = Vom S b W 


(36) 


An effective angle -of-attack is obtained for each sectional cut. The angle -of -attack 
imposed on the computational surface becomes 


<x c ( x ) = 


S b (X ) a 
S c (X) 


geom 


(37) 


By substituting three-point extrapolated differences for the velocities in the body 
boundary condition (Equation 32) and adding the thickness correction (35) and the 
angle -of-attack correction (37), the following expression for tody boundary points is 
obtained 



(38) 
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(38.1) 
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Solution Process 


Conventional finite difference relaxation schemes employ a procedure which is 
known as "Grid Halving". A solution is obtained using a very crude mesh. The grid 
is halved (mesh points are twice as dense) and the crude mesh potentials provide an 
initial flow field for what is now a medium mesh. The medium mesh solution is used 
for a mesh system which is halved once again for the final fine grid solution. Conver- 
gence is accelerated because flow disturbances naturally propagate faster in a crude 
grid than in a fine grid. 

The present method does not employ grid halving, however, an initial solution is 
obtained using only the global crude grid system. This quickly establishes the general 
character of the flow field and provides a starting flow for the detailed embedded grid 
systems. Final results are obtained by a simultaneous solution in all grid systems, 
both crude and fine. It is required then, that the different grid systems interact. This 
interaction is made possible by mesh overlap regions in which the flow' is computed 
twice during each cycle. 

Successive Line Over -Relaxation (SLOR) is used, vertical columns of mesh points 
are relaxed starting at the first point upstream and ending at the last point downstream 
in each grid system. The solution process for isolated bodies and wings will be 
described first followed by the solution process for wing -body combinations. 


Isolated Body Solution Process . 

Phase 1: The body is represented by crude mesh points and an initial solution is 

obtained in the global crude grid. Typically, 50 iterations are sufficient. 
Crude grid potentials (<f>) in the crude/fine grid overlap region (see Figure 7) 
provide a starting flow for the embedded fine grid. Linear interpolation 
is employed. 

Phase 2 (Step No. 1): The embedded body grid is swept holding fine grid perimeter 

points fixed forming an outer boundary. Body boundary points are computed 
and they form an inner boundary (see Figure 12). Body boundary potential 
values are computed using Equation (38). 
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(Step No. 2): The crude grid points which make up the inner boundary of the 
overlap region are updated based on fine grid potentials obtained in step 
No. 1. The global crude grid is swept holding these potential values fixed 
as an inner boundary. Infinity boundary conditions (Figure 1) at the limits 
of the computational space form an outer boundary. Crude grid potentials 
are now used to update fine grid perimeter points. Step No. 1 is repeated 

Phase 2 proceeds by repeating Steps 1 and 2 until both grid systems are con- 
verged. Typically, 50 cycles are sufficient. Body pressure coefficients are computed 
using body boundary potentials in the fine grid system. 


• • • BOUNDARY POTENTIALS FIXED 

/ REGION RELAXED 




STEP #1 

FINE GRID SOLUTION 


STEP #2 

CRUDE GRID SOLUTION 


Figure 12. Phase 2 Isolated Body Crude/Fine Grid Solution Process 
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Isolated Wing Solution Process. The isolated wing solution process is similar to 
that of the body. However, the entire fine embedded grid region is now the overlap 
region. In addition, side boundaries at the wing root and tip must be computed. 

Phase 1: The wing is represented by crude mesh points which fall randomly on the 

planform (see Figure 2). An initial solution is obtained using the crude 
mesh, typically 100 iterations are sufficient*. Using linear interpolation, 
crude grid potentials in the crude/fine overlap region (see Figures 3 and 4) 
provide a starting flow field for the fine embedded grid. 

Phase 2 (Step No. 1): The embedded wing grid is swept holding fine grid perimeter 

points fixed as an outer boundary (see Figure 13). Conventional Neumann 
boundary conditions ( IP ) are imposed at fine grid section boundary points 
forming an inner boundary. 



*See recommendations for usage section (pg. 148) for additional comments on 
iteration count and convergence. 
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(Step No. 2): The crude grid section boundary points are computed using the 
potentials at the fine grid section to undary points (linear interpolation). 

These crude potentials (v 3 ) are held fixed for the global crude grid sweep 
forming an array of dirichlet (<P) inner boundary conditions. Once again, 
infinity boundary conditions at the limits of the computational space form 
the outer boundary. At the end of the crude grid sweep, crude grid 
potentials are used to update the fine grid perimeter points. Step No. 1 
is repeated. 

Phase 2 proceeds by repeating Steps 1 and 2 until both grid systems are converged. 
Typically, 80 cycles are sufficient. Wing pressure coefficients are computed using wing 
boundary potentials in the fine grid system. 

Wing-Body Solution Process. The wing-tody solution process combines the 
procedures which are described for isolated bodies and isolated wings. In addition, 
however, a scheme for interacting two fine grid systems (embedded wing fine and body 
fine) is necessary. Figure 14 is a schematic illustrating the overlap of the three 
different grid systems involved. 

Phase 1: An initial solution is obtained with the wing -tody configuration represented 

in the crude mesh. Typically, 100 iterations are sufficient. Crude grid 
potentials are used to fill both wing and tody embedded grid systems. 

Phase 2: (Step No. 1) As was the case for isolated wings, the wing fine grid is swept 

using conventional Neumann boundary conditions at fine grid section bound- 
ary points. 

(Step No. 2) As was the case for isolated todies , the tody fine grid is 
swept with perimeter points and tody boundary points fixed. An additional 
surface representing the wing and wake is also held fixed (see Figure 15). 
Flow field potentials representing the wing/wake surface are obtained 
from the fine embedded wing grid. 

(Step No. 3) The global crude grid is swept with interior boundary points 
provided by the fine wing and tody grid solutions. At the end of the crude 
grid sweep, crude grid potentials are used to update fine tody grid perim- 
eter points and fine wing grid perimeter points outside the fine tody grid 
region (recall Figure 14). Within this region, fine tody grid potentials are 
used to update fine wing grid perimeter points. Step No. 1 is repeated. 
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Phase 2 proceeds by repeating Steps No. 1,2 and 3 until all grid systems are 
satisfactorily converged. Typically, 80 cycles are sufficient. 

The development of the solution in all cases is monitored by calculating the 
maximum update to the flow field potential ( <p~ -<p ) and its position in the three- 
dimensional flow field. Typically, for most configurations at flow conditions which 

are of interest, this maximum A v 3 value will start near 0. 5 x 10 1 and end near 
-4 

0.5 x 10 for the final solution. Large regions of supersonic flow will slow r this con- 
verging process since differ*ence approximations at supersonic points are only first 
order accurate w^hile those at subsonic points are second order accurate. Most impor- 
tant, however, the position of this maximum flow field correction will be useful in pin- 
pointing the problem area when convergence is hindered. Problem areas may develop 
when the flow conditions are extreme or when the geometric representation is in error. 


Body /Fuselage Geometry Model 

The present method provides a means for analyzing truly arbitrary body shapes 
at transonic speeds. Thus, it is occasionally required that complex three-dimensional 
geometries be input, processed and converted into a suitable array of boundary con- 
ditions for the analysis scheme. Since the input or modeling of complex body shapes 
is extremely error prone and certain applications might not warrant this level of effort, 
several options have been provided for describing body or fuselage shapes. 


Option No. 1 


Option No. 2 


Option No. 3 


Cylinder - A cylinder extending from upstream infinity to 
downstream infinity approximates the body shape. Only a 
radius need be input. This simple approach may be 
suitable for certain wing -fuselage applications (no body 
force or moment computed). 

Axisymmetric Body - Only an array of X and R coordinates 
need be input. This option is particularly advantageous when 
analyzing simple research test configurations. 

Arbitrary Fuselage - This option is necessary in most 
realistic aircraft applications , particularly when flow' details 
about fuselage contours (canopies, blisters and fairings) are 
required. 
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This section describes the mathematical modeling system which is employed for 
processing arbitrary fuselage shapes (Option No. 3). The system has been named 
"Quick-Geometry" by its developers, Vachris and Yaeger * 25 \ A detailed User's 
Guide for the Quick-Geometry System can be found in the Appendix of ref 26. This 
system was originally developed for the geometric modeling of wing -body shapes. 

Since only fuselage shapes are of concern here, many of the more sophisticated 
options including fillets and patches will not be described in the paragraphs which 
follow. In addition, if ref 25 and 26 are being used to augment the modeling descrip- 
tion provided herein, it should be noted that the input format has been modified to be 
more consistent with that of the basic transonic wing-body code. 

The geometry package requires that certain body lines and cross-section lines 
be defined. The body lines and cross-section lines may be likened to the stringers 
and bulkheads, respectively, used in fuselage construction. These line models are 
defined by a combination of simple curves (i.e. , lines, ellipses, cubics). They are 
taken together to provide a continuous analytical model of the surface geometry. 

Slopes and normals are developed analytically. Either discontinuous intersections or 
smooth fairings can be modeled and enforced. 

Two different coordinate systems are employed. Geometry definition is per- 
formed in a Cartesian coordinate system (x, y, z), while interrogation of the model for 
body boundary conditions is performed in cylindrical coordinates (X , R, © ). This 
results in the use of a plane of symmetry map axis , the height of which usually 
corresponds to the position of the max -half -breadth line. It is required that the 
configuration radius at any cross-sectional cut be a single valued function of the angle 
0 . These definition lines and coordinate systems are illustrated in Figure 16. 

From this description, it should be apparent that a minimum of four body lines 
are required for the simplest fuselage. These are 1) top centerline, 2) bottom center- 
line, 3) max-half -breadth line, and 4) the map axis. Each body line must be defined by 
both its Y and Z values over the full range of X (between fuselage nose and tail). 
Similarly, a minimum of two cross-section line segments are required for each 
different cross-section line model. These are 1) body upper, and 2) body lower. 

Both body lines and cross-section lines are specified by defining key arc or 
segment shapes and their accompanying limiters. The segment shape boundary 
conditions that are used to determine the coefficients of the shape equation are the 
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(1) origin point, (2) termination point, and (3) slope control point. The slope control 
point lies at the intersection of the line which is tangent to the segment shape at the 

origin point and the line which is tangent to the segment shape at the termination 
point (see Figure 17). The slope control point is a very convenient way of specify- 
ing slope conditions. In particular, it allows for the simultaneous specification of 
slope conditions at both ends of the segment. 



u 

Figure 16. Geometry Model Lines and Coordinate Systems 


36 



The arc shapes used for defining a cross-section line model are listed in Table 1. 


TABLE 1. Cross -Section Arc Shapes 


Shape 

Keyword 

Equation 

Line 

LINE 

Ay + B z + C = 0 

Ellipse 



(Concave to origin) 

ELLI 

2 2 
(y - y 0 ) + (z - z 0 ) =o 

A 2 B 2 

Ellipse 



(Convex to origin) 

ELLO 

Same as ELLI 


They are input in an order which starts at the body bottom centerline and proceeds to 
the body top center line. The segment shapes used for defining a body line model are 
listed in Table 2. 


TABLE 2. Body Line Segment Shapes 


Shape 

Keyword 

Equation 

Line 

LINE 

Ax + By = 0 

X-PARA BO LA 

XPAR 

Ax + By + y^ = 0 

Y-PARABOLA 

YPAR 

Ax + By + x 2 = 0 

X-ELLIPSE 

ELLX 

Ax + By + Cx 2 + y 2 = 0 

Y-ELLIPSE 

ELLY 

Ax + By + Cy^ + x^ = 0 

CUBIC 

CUBI 

Ax + By + Cx 2 + x^ = 0 


Cross-section arcs are input in their order of appearance. However, body line 
segments are defined along with an index which establishes their order in the x-direc- 
tion. In addition, body lines may be aliased to other body lines to avoid duplicate 
definitions . 

It should be noted that cross-sections are defined only in terms of named 
component arcs (arc shape table) and named control points. On the other hand, 
body lines are defined mathematically by coordinates over the length of the configura- 
tion for which they are required. At a given x-station, the body lines are interrogated 
to give the key control points required to construct the cross-sectional arcs. 
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Figure 17 is a schematic illustrating the component build-up of a particular 
body line and cross-section line model. Naturally, LINE segments do not require a 
slope control point. In this case, the portion of the body top center line illustrated 
requires four body line segments and the cross-section is constructed with two arcs 
(two is the minimum number allowed). 

Figure 18 illustrates an array of cross-section lines and the body lines for a 
typical transport fuselage shape. Five different body lines can be identified. These 
body lines form the limiters for the three distinct cross-section models. This geometry 
package has proven to be particularly useful when coupled with the mesh embedding 
scheme since the fine grid structure may be placed for analysis purposes about any 
region of interest in addition to totally encapsulating the body shape. Figure 19 
illustrates the surface geometry produced when the mathematical model is interrogated 


COMPONENTS OF BODY (TOP CENTER) LINE 
MODEL. SEE FIG. 16. 


CANOPY SLOPE 
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for a fine embedded mesh system. Boundaries of this region were outlined in Figure 
18. Note that the mathematical model does not change in this case, only the region 
over which the model is interrogated for boundary conditions and the surface shape. 

The input data set for this geometry can be found in the Appendix (page 109 and 
computed results can be seen in the Results Section (page 59). 

Figure 20 illustrates an array of cross-section lines and the body lines for an 
area-ruled transport fuselage shape. Here, only four body line models and one cross- 
section model were required. The cross-sectional shape has been varied from circular 
to flat bottomed simply by changing the position of the max -half -breadth line. 

The geometry modeling procedures described in the preceeding paragraphs can 
best be interpreted by studying a sample fuselage shape. Table 3 lists the geometric 
model for the transport nose -winds hie Id -canopy shape illustrated in Figure 19. Lines 
C through S make up the cross-section line model input while lines T through FFF 
describe the body line models. Note that there are three distinct cross-section line 
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Figure 18. Crude Grid Interrogation of Complete Geometry Model 
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TABLE 3. Gulfstream Fuselage Quick-Geometry Model 
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models (lines C-E, F-I, and J-N). Line group P-S provides the region over which 
each cross-section line model is to be applied. Most important however, note that 
line S indicates that the same cross-section model is to be used for both the nose and 
fuselage tody. 

Lines T through FFF indicate that ten body lines are required for this fuselage 
shape. Each is described in both y and z values over the range of x for which they are 
required. This leads to twenty different models along with the y center line used for 
aliasing. Ten are defined by inputting section shape and limiters (lines U-UU) and 
eleven are aliased to them (lines VV-FFF). 


As one might expect, it is very easy to make errors when defining complex 
three-dimensional shapes. Various diagnostic checks have been incorporated into 
the method to simplify this process. First, geometry verification plots are available 
to illustrate the defined surfaces using cross-section lines. Second, a diagnostic print- 
out is available for checking both the coefficients of the shape equations and the proper 
range of definition for tody segments. Finally, the sample cases provided in the 
Appendix should provide a foundation of experience upon which a user can build. 

Wing Viscous Effects 


Wing viscous effects (C $ *) are computed in the present method by coupling 
a modified Bradshaw boundary layer computation with the aforementioned finite 
difference potential flow scheme. The boundary layer method provides details of the 
thin viscous layer close to the surface given the wing pressure distributions from the 
global inviscid calculation. Several modes of operation are provided for tailoring the 
computation for various applications and minimization of computing requirements. 


Mode No. 1 Wing viscous effects including skin friction drag component 

(Cj)j) and boundary layer displacement thickness ( J *) are not 
computed. Inviscid solution only. Minimum computer resources 
required. 

Mode No, 2 Wing viscous effects are computed at the end of the inviscid 
solution providing an approximation to and 5 *. Permits 
subtraction of B. L. 5 * for design work done inviscidly to get 
physical or constructed wing shape (approximately 7 % computing 
time increase over Mode No. 1). 
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Mode No. 3 Coupled inviscid/viscous interaction solution. Iterative scheme 
for simultaneous solution of boundary layer and potential flow 
(approximately 30% computing time increase over Mode No. 1). 

Examples of viscous flow phenomena occurring on wings have been sketched in 
Figure 21. In particular , an aft-loaded or supercritical wing section shape has been 
illustrated. Details of the shock wave /boundary layer interaction are not computed 
by the present method. This effect includes both a local thickening of the boundary 
layer and a weakening of the shock w’ave strength. Instead, non -conservative finite 
differences of the flow' equation have been selected which somewhat fortuitously 
approximate the shock weakening phenomena. The global boundary layer displacement 
thickness effects are predicted by the present method. For aft-loaded wing sections 
or wings with control surface deflections , this effect can have a dramatic effect on 
wing loading and shock wave positions. As illustrated in Figure 21, an asymmetric 
boundary layer build-up results in a net decambering effect. For Mode 3 operation, 
this effect is incorporated into the solution automatically. 
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(2 7 ) 

The modified chord technique of Nash & Tseng is used in the present method. 
Developed by Mason ^ 28 ', this particular technique permits a two-dimensional boundary 
layer method to be extended to the three-dimensional case provided that the flow does 
not deviate far from the infinite sheared wing type. The coordinate system for this 
sheared wing approximation can be seen in Figure 22. The wing sw T eep angle is taken 
to be that of the mid-chord span line. As a result, the local boundary layer computa- 
tions vary with the wing planform shape, see Figure 23. This quasi-three-dimensional 
approach has been selected because it requires about an order of magnitude less 
computing time as compared to a full three-dimensional boundary layer calculation. 
This is an important consideration when coupling with three-dimensional transonic 
flow methods is being investigated. 
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MID-CHORD 
SPAN LINE 



(29) 

The Bradshaw turbulent boundary layer method ' ' provides the foundation for 

this scheme with the laminar boundary layer predicted by Thwaites method modified 

for compressible flow. Several corrections are required for simulation of three- 

( 27 ) 

dimensional effects. The modifications of Nash and Tseng for incompressible 
flow about an infinite yawed wing are employed. In particular, it was noted that, 

A) pressure gradient effects on boundary layer development were properly accounted 
for by a 2-D calculation performed in the chordwise x" direction and similarly, B) shear 
stress direction and magnitude could be represented by a 2-D calculation in the 
streamline direction (see Figure 22). Nash has shown that the equation for the 
chordwise shear stress component is similar to the conventional 2-D equation. 

T x = cos [tan-1 ^ e ] r = Xt (39) 

Here it is assumed that X varies slowly in the 3c-direction compared to the shear 
stress change. 
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The computation requires an effective Reynolds number which is obtained by 
requiring that the shear component in the x-direction bear the same relationship to 
an effective Reynolds number that the actual shear stress has to the specified Reynolds 
number. This is accomplished by using the Reynolds number based on the momentum 
thickness at transition in the Karman-Schoenherr formula. The transition point has 
been fixed at the 5 % chord position in the coding of the method. Of course, this may 
be modified or varied by changing the coding. However, this variation is not provided 
as an input option. 

Finally, a method is required to compute the solution in the spanwise direction. 

It is assumed that r~can be related to 0 by a flat plate skin-friction formula. In 

Z 

this case, the Squire-Young formula is used. 

r w^ l 

p w ^ d'y 9 

00 [5.89 log-, n (4.075 Re — )] 2 

iU c (40) 

NOTE : Re based on W«> and effective chordwise viscosity of Nash. 

The boundary layer computation will stop when flow separation is predicted 
( T = 0). At this point, the slope of the displacement thickness boundary is 
extrapolated to the trailing edge. This will permit the calculation to proceed, however, 
it should be noted that deviations between computed results and the real flow will 
result if separation points before the 96 to 97% chord position are encountered. 

The interaction solution (coupled viscous/inviscid) requires that the boundary 
layer computation be performed during the solution development. The computed 
boundary layer displacement thickness slopes are added to the wing geometric surface 
slopes to provide an equivalent inviscid wing shape for analysis by the inviscid 
scheme. Note that in the present small -disturbance method, only the boundary con- 
dition slopes are modified. The planar wing boundary surface does not change. 

Many wing shapes which are of practical interest have section shapes with cove 
regions. This includes the supercritical type wing section shapes and occasionally 
conventional sections when control surfaces are deflected. Of course, it is possible 
for cove type shapes to exist on the v'ing upper or lower surface. The existence of 
this type of geometry coupled with extreme flow conditions can become a problem 
during the solution process if cove type flow' separations are encountered. Typically, 
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the boundary layer becomes thicker in the cove region and thins out toward the 
trailing edge (see Figure 21). If separation is encountered, the boundary layer com- 
putation will stop and the character of the displacement thickness effect will be lost. 

In this case , the section lift may be significantly in error which will lead to a poor 
numerical simulation. To relieve this problem, an empirical relation developed by 
Bavitz is employed to provide an approximation to the cove displacement thick- 
ness. This will make it possible for the computation to continue. 

Bavitz's empirical relation is based on an extensive correlation study of 
supercritical type airfoils. It involves the pressure coefficients and displacement 
thickness at four points along the section surface. These points have been illustrated 
in Figure 24. The first point is positioned 109? chord upstream of the separation point. 
The second point is 89? chord before separation, the third point is midway between the 
separation point and the trailing edge, and the final point lies at the trailing edge. 


The boundary layer displacement thickness at stations 3 and 4 are computed using 
using the following equations: 
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(41) 



(42) 


A third degree polynomial is fit to the four 8 * values to provide the boundary layer 
shape between the separation point and the trailing edge. The pressure difference 
coefficient in Equation (41) was set at 0.033 for computations involving two-dimensional 
airfoils. Calculations performed on several supercritical wings, however, indicated 
that a lower value would be better for the three-dimensional case. As a result, the 
0.019 coefficient is used in the present method. 

It should be noted that this empirical relation only provides an approximation 
to the cove region boundary layer when shallow separations are encountered. In many 
cases, this approximation may be sufficient to establish the proper wing lift level and 
the final solution in the cove region may be fully attached (empirical correction not 
used). If this is the case, a valid solution has probably been obtained. If on the other 
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hand, the final solution still indicates cove separation is present, it is likely that the 
physical flow is separated and agreement between theory and experiment may be poor. 

In the present method, the empirical correction is employed if separation is 
encountered and the trailing edge angle is less than -4° or greater than +4°. For 
trailing edge angles within these boundaries, the conventional procedure is used beyond 
the separation point. In other words, the sum of the o * and local wing slopes at 
the separation point is enforced as the wing boundary condition between the separation 
point and the wing trailing edge. Viscous effects are incorporated by using an under- 
relaxation factor of 0.6. This factor should be reduced if boundary layer instabilities 
are encountered. 
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Pressure, Force, and Moment Coefficients 

Force (except wing friction) and moment coefficients are obtained by integrating 
computed pressure coefficients. Wing pressure coefficients are computed by using 
central difference operators for velocity components (Equations 17a, b) and then 
pressure coefficients from Equation 2. With 100 grid points between each section 
leading and trailing edge, this provides pressure output every l c /< chord beginning 
at 0.005c. 


Body pressure coefficients are computed in a similar manner except that a 
correction must be applied in regions where the computational body surface and the 
true body surface are not aligned. The velocity or pressure distribution corrections 
are similar to those developed for the body boundary conditions (Equation 38) . Both 
apply only in regions where there is a displacement between the computational body 
surface and the actual body surface. Velocity components on the computational sur- 
face are obtained with central difference operators at body boundary points. The 
velocity on the true body shape is then approximated by 


where 


u b = ^x 



+ a7 ( ^ 5r 


(43-A) 


5r= R b -R c 


(43-B) 


Since the flow is continuous in this region, (Equation 43-A) may take the form 






+ JX ( $r> 5r 
c c 


+ "^n Sr 
c c 


Using (Equation 35) the body surface velocity becomes 


(44-A) 


(44-B) 


u b u c + ( R b R c> dX ^n 

c 


(45) 
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Body surface pressures are now computed using a central difference operator 
(Equation 17a) for the surface velocity (Equation 45) which is substituted into the 
pressure equation (Equation 3). 


Configuration force , moment , and loading coefficients 1 ’equire reference 
areas and lengths as a basis. These reference values are taken from the input 
geometric wing planform shape. A typical wing planform has been sketched in 
Figure 25. Note that the wing is specified to the plane of symmetry even for the 
wing -body case. The following parameters are computed using the defined wing 
planform . 
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( 48 ) 
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Wing section coefficients are obtained by integrating the upper and lower 
surface pressure coefficients along the chord. 
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In addition, a sectional moment about the axial reference position is computed. 


Cm cj M.A.C. C , J (C P U C P p )< X loc X REF )dx 

loc -b u * ( 52) 

The wing lift coefficient is obtained by integrating the spanwise load coefficients. 

c -2 C W I C|ocC|? “ ) H 

L WING b J o 1 C, / dy (53) 

Note that for wing -body combinations, all integrals in the spanwise direction start 
at the wing -body juncture instead of the symmetry plane which is the case for an iso- 
lated wing. In addition, the computed spanwise integrals are scaled by S /S for 
the wing -body case. 

The wing moment coefficient is computed in a similar manner. 


2 f h/ 2 ( C loc C m w ^ 
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Wing skin friction drag is computed by integrating the shear stress coefficients 
computed by the boundary layer analysis. At each wing station, a local section skin 
friction coefficient is computed. 




-f 


[ T x cos \oc + T z Sin A locl dx 


The x and z sheared wing coordinate system was illustrated in Figure 22. Total \\ing 
skin friction drag is then obtained by integrating the local friction coefficients aLong 
the span. 

C =1 
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The wing pressure drag coefficient: 
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is added to the skin-friction drag coefficient to obtain the total wing drag coefficient. 

C °WING " Cf WING + C ° P(W) (58 ) 

Body coefficients are obtained by integrating surface pressure distributions and 
adding a skin friction component based on total body wetted area. A body cross- 
sectional force coefficient distribution is computed first. 
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The longitudinal coefficients are then integrated along the length of the body. 
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The body moment coefficient (62) is scaled by L/M.A.C. for the wing-body case. 

The reference position about which moments are computed is taken to be the 
center of the body for isolated bodies and the input moment center for the wing -body 
case. In addition, the integrated coefficients are based on reference areas which are 
a function of the case. For isolated bodies , the maximum cross-sectional area 
becomes the reference area while for wing -body configurations, the wing plan form 
reference area is used. For isolated bodies, equations (61-63) are scaled by 
Sr / s b aq - and for wing -body combinations, the equations are scaled by S Bp /S w 

The body skin friction coefficient is computed using the Prandtl-Schlichting 
formula 

r _ = 0.455 

F BODY _ [lo g (Reg)] 2-58 

(64) 

corrected for compressibility effects. Here, Re^ is the Reynolds number based on 
body length. Total body drag is then computed. 

d BODY F BODY D P(b) (65) 


53 



Total configuration force and moment coefficients are obtained by adding the wing and 
body components. 


° L Cl WING + Cl BODY 
C M C M W j NG + C M B q DY 

Cp) ~ C + Cp. 

u d wing d body 

Configuration lift induced drag (C Dl ) is computed by first determining the span- 
wise loading efficiency factor (E). This is accomplished by means of fast Fourier 
analysis applied to the configuration span load distribution. The lift induced drag 
component is then computed using the following equation. 

c l 

C Di jr /RE 

(67) 

The wave drag is then computed by subtracting the lift induced drag from the total 
integrated pressure drag given by (Equations 57 and 63). 


(66a) 

(66b) 

(66c) 


C °WAVE CD P(co) + C Dp (b ) % (g8) 

It is important to note that the body, wing, and wing -body force and moment 
coefficients are based on reference areas and lengths taken from the input geometry. 
These values are listed as part of the program output. Sometimes it is required 
that configuration force amd moment coefficients be based on reference values 
other than those provided by the physical wing and body shape. If this is the case, and 
case 8 is a good example, it will be necessary to scale the computed coefficients by 
the ratio of the two reference values. If A is the reference length or area, the actual 
coefficients are computed. 
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COMPARISONS AND TYPICAL RESULTS 


The computational method is capable of treating a variety of geometric shapes. 
Several options are available for modeling and analyzing a particular configuration. 

If a single sample case were chosen as a demonstrator, many features of the method 
would not be illustrated. Since transonic experimentation is often impaired by wind 
tunnel wall interference effects, it would also be bad practice to rely on a single case 
to evaluate a numerical method. For this reason, eight different sample cases have 
been selected. These include three isolated body or fuselage shapes, two isolated 
wings, and three wing -body combinations. These cases cover both sharp and blunt nose 
bodies as well as conventional and supercritical type wing sections over a range of flow 
conditions . 

The comparisons are of particular interest because the present method is based 
on certain small-disturbance assumptions. It is important then to study these applica- 
tions if the method is to be effectively applied to realistic- aircraft configurations. The 
input data set for each of these cases has been included in the appendix of this report. 
Figure 26 defines the various symbols used for the correlations throughout this 
section. 

Isolated Bodies 

A sharp -nose body of revolution was investigated by Swihart and Whitcomb (31). 

The body had a fineness ratio of 12 and was sting-mounted. This data was selected 
because it provided isolated body experimental results for non-zero angles -of -attack. 
Correlations for the zero angle -of -attack case at M = 0.99 can be seen in Figure 27. 

A computation was also made for the 8.4° angle -of -attack case at M = 0.99. Correlation 
with experimental data for the top and the bottom centerlines can be seen in Figure 28. 
Discrepancies on the upper surface near the body -sting juncture can probably be attributed 
to boundary layer thickening. 

Experimental data for a blunt -nose body of revolution was recently obtained by 
Couch and Brooks This body had a fineness ratio of approximately 9. The cross- 

sectional area distribution is typical of supercritical type bodies. The pressure 
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Figure 26. Results Section Symbol Definition 
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Figure 27. Sharp Nose Body Pressure Distribution Correlation 
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Figure 28. Sharp Nose Body Pressure Distribution Correlation 


DEGREES 



distributions are characterized by a rapid expansion around the nose followed by a flat 
distribution along the length of the body. Correlations for the case M = 0.980 and 
M = 0.991 can be seen in Figure 29. 

These simple research models for which comparisons were just presented pro- 
vide a means of checking the behavior of the present modified small disturbance 
method when analyzing both blunt and sharp-nose bodies. However, these shapes are 
rarely of interest to the aircraft designer. To demonstrate the full computational 
capability of the present method, calculations have been made for a realistic transport 
fuselage. The geometry for this case can be seen in Figures 18 and 19. In this par- 
ticular application, the fine mesh is positioned for detailed analysis of the nose-wind- 
shield-canopy region of the fuselage. 

Wind tunnel pressure distributions were measured to provide an understanding 
of shock-wave /boundary layer interactions over the canopy which caused unacceptable 
levels of noise in the cockpit. The sting-mounted isolated fuselage model was 
instrumented to obtain detailed pressure data in the windshield-canopy region ( 33 ). 

The objective of follow -up experiments or the use of the analysis method would be to 
re-contour the fuselage shape to reduce the strength of the canopy shock wave. Corre- 
lations with experimental data for the fuselage top centerline at the transonic cruise 
condition can be seen in Figure 30. The computed decay of the canopy shock wave can 
be seen in Figure 31. 
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Figure 29. Blunt Nose Body Pressure Distribution Correlation 
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Figure 30. Grumman Gulf stream II Fuselage Pressure Distribution Correlation 
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Isolated Wings 

Two different wings have been selected for comparisons. The first is a low' 
aspect ratio wing with a symmetric wing section and the second is a high aspect ratio 
wing with a conventionally cambered NACA 64^-212 wing section. Viscous effects on 
the symmetric wing are negligible and no viscous effects w'ere computed (VISMOD - 1). 
Viscous effects did not have a large impact on computed pressures for the cambered 
wing but they did significantly change the computed forces and moments. The cambered 
wing was analyzed using the viscous /inviscid solution process (VISMOD = 3). 

The ONERA M6 wing (34) provides a very simple case for the computational 
method because wing planform sweep and taper are mild. In addition, viscous effects 
are negligible making this a very popular case for published correlation studies. The 
ONERA wing is interesting, however, because the data exhibits a double shock w F ave 
pattern. Superimposed computed pressure distributions for the w T ing upper and lower 
surface can be seen in Figure 32. A weak highly sw'ept supersonic to supersonic shock 
wave forms at the wing leading edge. As the wing tip is approached, this w^eak shock 
w'ave coalesces with the primary supersonic to subsonic shock wave behind it. Published 
reports indicate that methods with poor resolution or flow' governing equations without 
additional cross ’-flow 1 terms cannot resolve the weak highly sw^ept type of shock wave. 
Correlations with experimental data at five span stations can be seen in Figure 33. 
Spanwise distributions of lift, drag and pitching moment can be seen in Figure 34. 

Figure 35 shows the computed spanload and drag component breakdown. 

A cambered high aspect ratio wing (^ w 7 as analyzed using the present method. 

This case is thought to be of particular interest because the geometry and lift condition 
are typical of conventional transport wrings now in existence (super-critical type wings 
are also included - case number 8). Figure 36 illustrates the wing shock wave pattern. 
Correlations at five spanwise wing stations can be seen in Figure 37. The computed 
spanload along with span efficiency factor and drag component breakdown can be seen 
in Figure 38. Force and moment comparisons can be seen in Figure 39. 
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Figure 31. Grumman Gulf stream II Fuselage-Decay of Canopy Shockwave, M = 0.800, a * 3.10 c 
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Figure 32. Superimposed Computed Pressure Distributions for isolated Onera M6 Wing 
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Figure 33. Onera M6 Wing Pressure Distribution Correlation (Sheet 1 of 2) 
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Figure 33. Onera M6 Wing Pressure Distribution Correlation (Sheet 2 of 2) 
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Figure 35. Onera M6 Wing Computed Spanload 
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Figure 36. NACA RM A9K01 Isolated Wing Superimposed Computed Pressure Distributions 
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Figure 37. NACA RM A9K01 Wing Pressure Distribution Correlation (Sheet 1 of 2) 
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Figure 37. NACA RM A9K01 Wing Pressure Distribution Correlation (Sheet 2 of 2) 
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Figure 38 NACA RM A9K01 Isolated Wing Spanload Correlation 







Wing Body Configurations 

Three different wing -body configurations have been selected for correlation 
studies. They are 1) a simple wing-body research model with a planar wing, 2) a 
conventional wing transport aircraft, and 3) a supercritical w'ing area-ruled fuselage 
transport configuration with wing control surface deflection effects. The final con- 
figuration with its highly aft -cambered wing section illustrates both the use of the 
viscous/inviscid interaction mode of operation and the arbitrary body modeling scheme. 

A simple wing-body research model < 36) is illustrated in Figure 40. The center 
mounted wing has 45° of sweep with a NACA 65A006 wing section. Superimposed 
computed pressure distributions for the wing can be seen in Figure 41. Correlations 
with experimental wing pressures at five spanwise stations can be seen in Figure 42. 
Body pressure correlations can be found in Figures 43 and 44. Boundary layer dis- 
placement thickness effects were not computed for the symmetric wing in this case. 

Force and moment comparisons are made in Figures 45 through 49. The configuration 
drag source breakdown can be seen in Figure 48. 

The Boeing KC-135 transport has been illustrated in Figure 50. This is a 
good example of a conventionally winged transport. Comparisons have been made for 
the configuration without engine pods or pylons. An initial computation indicated that 
this wind-tunnel model may be experiencing a wing twist caused by wing loading. As a 
result, the geometry analysed includes 1° of negative section incidence at the wing tip 
blending to 0° at the wing -body juncture. Superimposed computed wing pressures can 
be seen in Figure 51 for the cruise design point. Correlations at three wing span 
stations can be seen in Figure 52. 

The final case is a supercritical wing/area ruled fuselage transport config- 
uration (See Figure 53). The solution diverged if the actual glove leading edge 
sweep of 72° was modeled. For this reason, a reduced sweep of 57° was employed for 
these comparisons. It is expected that the wing root leading edge pressure peak 
discrepancy (Figures 55 and 56) is a result of this approximation. In addition, the wing- 
body computed moment coefficient will have an erroneous nose down component (Figures 
58 and 59). The fuselage model for this calculation was shown in Figure 20. This 
model is interesting because it illustrates the aerodynamic effects that a wing control 
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surface deflection has on a supercritical wing when the flow is transonic and shock 
waves are present. For these comparisons , data at a -15° deflection of the second 
control surface was used. The control' surface can be seen in Figure 53 and on the 
wing section pressure distribution plots. Superimposed computed pressure distri- 
butions at the cruise design point can be seen in Figure 54 for the basic wing shape. 
Wing pressure distribution correlations at five spanwise stations can be seen in Figure 
55. A similar set of comparisons for the wing with the control surface deflected up- 
ward can be seen in Figures 56 and 57. Force and moment comparisons are made in 

Figures 58 and 59 while the spanload and drag source breakdown can be seen in 
Figure 60. 
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Figure 40. NACA RM L51F07 Simple Wing-Body Configuration 
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NACA RM L51F07 AX I -STM B0DT/4S DEG SWEPT WING 
WING PRESSURE DISTRIBUTION - UPPER SURFACE 
MACH = 0.930 ALPHA = >4.00 RE = 2.03 M 

WING CL = 0 . 2S8 CM = -0.048 CO = 0.0230 


NACA RM L51F07 AXI-5TM BODt/45 DEG SWEPT WING 
WING PRESSURE DISTRIBUTION - LOWER SURFACE 
MACH = 0.930 ALPHA = 4.00 RE = 2.02 M 

WING CL =0.258 CM = -0.048 PH = 0.0230 


Figure 41. NACA RM L51F07 Wing-Body Superimposed Computed Wing Pressure Distributions 
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NflCfi RM L51F07 AXI-STM BODT/45 DEG SWEPT WING 
WING STATION 5 2T/B = 0.186 

MACH * 0.930 ALPHA =4.00 

SECTION CL = 0.293 CM = -0.051 CO = 0.0310 


NACA RM L51F07 AX 1 -STM BOOT/45 DEG SWEPT WING 

WING STATION 9 2T/B = 0.382 

MACH = 0.930 ALPHA =4.00 

SECTION CL = 0.345 CM = -0.039 CD = 0.025 


Figure 42. NACA RM L51F07 Wing Pressure Distribution Correlation (Sheet 1 of 2) 
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NflCR RM L51F07 RXI-5TM BOD 1745 DEG SWEPT WING 
WING 5TRTI0N 13 2T/B = 0.603 

MRCH = 0.930 RLPHfl = 4.00 

SECTION CL = 0.364 CM = *0.009 CD = 0.0199 


Vw‘ 

*0 


V = 0.80 



NRCR RM L51F07 RXF-5TM B00T/45 DEG SWEPT WING 
WING STATION 16 21 / B = 0.799 

MRCH = 0.930 PLPHR = 4.00 

SECTION CL - 0.334 CM = 0.021 CO = 0.0146 



NHCfl RM L5IF07 RX1-STM BODT/45 DEO SWEPT WING 
WING STATION 18 21/B = 0.956 

MRCH = 0.930 RLPHR = 4,00 

SECTION CL = 0.211 CM = 0.033 CO = 0.01 


Figure 42. NACA RM L51F07 Wing Pressure Distribution Correlation (Sheet 2 of 2 ) 
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Figure 43. NACA RM L51F07 Stations for Computed and Experimental Body Pressures 
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NfiCfi RM L51F07 fix I -STM B00Y/4S DEG SWEPT WING NflCH RM LSI -07 flXl-STM BGDT/45 OEG SWEPT WING 

BOOT PRESSURE DISTRIBUTION BODt PRESSURE DISTRIBUTION 

BODY STATION 1 BODY ANGLE = 90.00 DEGREES BOOT STATION 6 BODT ANGLE = 45.00 DEGREES 

MflCH = 0.930 ALPHA = 4.00 MACH = 0.93G ALPHA = 4.00 



NACA RM LS1F07 fix I -SIM B0DT/4S DEG SWEPT WING 

BODT PRESSURE DISTRIBUTION 

BODY STATION 9 BODt ANGLE = 17.35 DEGREES 

MACH = 0.930 ALPHA = 4.00 
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NACA RM L5IF07 AX I -STM BODT/45 DEG SWEPT WING 
BOD i PRESSURE DISTRIBUTION 

BODY STATION 20 BOOT ANGLE = -17.35 DEGREES 
MACH = 0.930 ALPHA =' 4.00 


Figure 44. NACA RM L51F07 Body Pressure Distribution Correlation (Sheet 1 of 2) 
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NACA RM L51F07 AX I -STM BODT/45 DEG SWEPT WING 

BOOT PRESSURE DISTRIBUTION 

BODY STATION 17 BODY ANGLE = -45.00 DEGREES 
MACH = 0.930 ALPHA = 4.00 
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NACA RM L51F07 fix I -STM BODY/45 DEG SWEPT WING 
BODY PRESSURE DISTRIBUTION 

BODY STATION 1 2 BODY ANGLE = -90.00 DEGREES 
MACH = 0.930 ALPHA = 4.00 


Figure 44. NACA RM L51F07 Body Pressure Distribution Correlation (Sheet 2 of 2) 
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NACA RM L51F07 AXI-STM BODY/45 DEG SWEPT WING 

SPAN EFFICIENCY E =0.927 

LIFT INDUCED DRAG =0.0072 

WAVE DRAG = 0.0094 

FRICTION DRAG =0.0135 

Figure 45. NACA RM L51F07 Wing-Body Spanload Correlation 


75 



E 




body longitudinal lord cl r/rmrx 



0- f MO o.'uoo o' 600 o'eoo Tooo 

BODY POSITION X/L 

NACA RM LS1F07 AXI-SYM B0DT/4S DEG SWEPT WING 

BOOT LONGITUDINAL LOADING 

BODY CL = 0.0319 

BODY CM - 0.1325 

BODY PRESSURE CO = 0.0000 

BODY FRICTION CD = 0.0072 


Figure 49. NACA RM L51F07 Body Longitudinal Load Correlation 
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BOEING KC-135 TRANSPORT 
WING PRESSURE DISTRIBUTION - UPPER SURFACE 
MACH = 0.780 ALPHA =3,00 RE = 5.25 M 

WING CL =0.330 CM = -0.010 CD = 0.0122 


POE INC KC-1^ TRANSPORT 
Wlf-iG PRESS -:.' DISTRIBUTION - LOWER SURFi 
MACH = C.7EC Alpha =2.00 PE = I 
WING CL = 0.330 CM = -0.010 CD = 


Figure 51. Boeing KC-135 Transport, Superimposed Computed Wing Pressure Distributions 
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wing S'r: k- e 21/8 = i .332 

MACH = 0,7 b: ALPHA = 2.00 

SECTION CL - 0 . 46c? C:i - - 0.030 CD - 0.0157 


BOEING KC-135 TRANSPORT 

WING STATION 16 2T/B = 0.799 

MACH = 0.780 ALPHA = 2.00 

SECTION CL = 0.601 CM = -0.0M4 CO = 





BOEING KC-135 TRANSPORT 

WING STATION 16 2Y/B = 0.956 

MACH = 0. 780 ALPHA = 2.00 

SECTION CL * 0.534 CM = -0.014- CD = O.OOMO 

Figure 52. Boeinq KC-135 Transport Wing Pressure Distribution Correlation 
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NASA TM X-3431 TRANSPORT CONFIG. 9 FEB 1979 NASA TM X-3431 TRANSPORT CONFIG. 9 FEB 1979 

WING STATION 5 2T/B = 0.186 WING STATION 7 2T/B = 0.282 

MACH = 0.900 ALPHA = 3.91 MACH = 0.900 ALPHA = 3.91 

SECTION CL = 0.412 CM = -0.068 CD = 0.0368 SECTION CL = 0.636 CM - -0.102 CD = 0.0327 



NASA TM X-3431 TRANSPORT CONFIG. 9 FEB 1979 NASA TM X-3431 TRANSPORT CONFIG. 9 FEB 1979 

WING STATION 9, 2T/B = 0.382 WING STATION 12 . 2T/B = 0.544 

MACH = 0.900 ALPHA = 3.91 MACH = 0.900 ALPHA = 3.91 

SECTION CL = 0.659 CM = -0.128 CD = 0.0501 SECTION CL = 0 . 740 CM = -0.173 CD = 0.0375 

Figure 55. NASA TM X-3431 Transport Configuration Wing Pressure Distribution 
Correlation for Basic Wing (Sheet 1 of 2) 
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NASA TM X- 34 3 1 TRANSPORT CONFIG. 9 FEB 1979 
WING STATION 15 2T/B = 0.730 

MACH = 0.900 ALPHA = 3.91 

SECTION CL = 0.703 CM = -0.145 CO = 0.0158 

Figure 55. NASA TM X-3431 Transport Configuration Wing Pressure Distribution 
Correlation for Basic Wing (Sheet 2 of 2) 



§J NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL OEF. 
J WING PRESSURE DISTRIBUTION - UPPER SURFACE 
MACH = 0.900 ALPHA =3,93 RE = 1.58 M 

WING CL =0.327 CM = -0.054 CD = 0.0249 



Figure 56. NASA TM X-3431 Transport Configuration Superimposed Computed Wing 
Pressures for Wing with Control Surface- Deflected 
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NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 
WING STATION 5 3T/B = 0.186 

MACH = 0.900 ALPHA = 3.93 

SECTION CL = 0.342 CM = -0.054 CO = 0.0291 



NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 
WING STATION 7 2T/B = 0.282 

MACH = 0.900 ALPHA = 3.93 

SECTION CL = 0.291 CM = 0.014 CO = 0.0152 
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NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 
WING STATION 9 2Y/B = 0.382 

MACH = 0.900 ALPHA = 3.93 

SECTION CL = 0.452 CM = -0.063 CO = 0.0401 
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NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 
WING STATION 12 2T/B = 0.544 

MACH = 0. 90 C ALPHA = 3.93 

SECTION CL = 0,603 CM = -0.150 CD = 0.0316 


Figure 57. NASA TM X-3431 Transport Configuration Wing Pressure Distribution 
Correlation for Wing with Control Surface Deflection (Sheet 1 of 2) 
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NASH TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 
WING STATION 15 2T/B = 0.730 

MACH = 0.900 ALPHA = 3.93 

SECTION CL = 0.644 CM = -0.165 CD = 0.0209 


Figure 57. NASA TM X-3431 Transport Configuration Wing Pressure Distribution 
Correlation for Wing with Control Surface Deflection (Sheet 2 of 2) 
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NASA TH X-3431 TRANSPORT CONFIG. 9 FEB 1979 

SPAN EFFICIENCY E =0.911 

LIFT INDUCED DRAG =0.0116 

WAVE ORAG = 0.0106 

FRICTION DRAG = 0.0155 


NASA TM X-3431 TRANSPORT CONFIG. WITH CONTROL DEF. 

SPAN EFFICIENCY E =0.722 

LIFT INDUCED DRAG =0.0092 

WAVE DRAG = 0.0084 

FRICTION DRAG - 0.0155 


Figure 60. NASA TM X-3431 Transport Configuration Effect of Wing Control Surface 
Deflection on Computed Wing Spanload 
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CONCLUDING REMARKS 


The present method provides certain advantages in both the design and analysis 
of wing-fuselage configurations at transonic speeds. These advantages are related to 
the high density resolution available and the grid system flexibility needed for treating 
complex three-dimensional shapes. However, to make these gains, some desirable 
features had to be sacrificed. In particular, the basic small disturbance-planar 
boundary condition approach exhibits certain limitations when compared to the more 
sophisticated full potential equation - exact boundary condition approach. This method 
can be used effectively if these limitations are understood. 

It was pointed out in the grid system section that wing section leading edges are 
positioned between fine grid points. It is not possible to have a boundary point exactly 
at the leading edge, a restriction inherent in the small disturbance formulation. As a 
result, details of the flow in the wing section nose region are somewhat compromised. 
Unfortunately, the drag computations involve a pressure integration that includes the 
nose region where for blunt sections, a suction force naturally develops. For this 
region, it is expected that the present method will be useful in predicting drag 
increments but caution should be exercised when using the absolute drag levels. 

It should be noted that for certain extreme geometric configurations and certain 
extreme flow conditions, computed small disturbance solutions may not correlate 
with the physical flow as well as a full potential equation solution. The Douglas wing- 
body configuration (Ref 39) is a good example (see Figure 61). This wing has 15% 
thick blunt wing sections. If the flow condition is mild, for example, below the drag 
divergence condition, comparisons with experimental data are good (see Figure 62 
and 63). However, if flow conditions are extreme, beyond what would be considered 
the design point or drag divergence conditions, comparisons are somewhat compromised* 
(see Figure 64). However, computations indicate that the present method performs 
well for the more common geometries and flow conditions occurring in most practical 
applications. In particular, existing transports analyzed at normal cruise conditions 
have yielded excellent correlation with experimental data. To date, these transports 
include the Boeing KC-135, Grumman Gulfstream III, and Lockheed C-5A and C-141. 


*See Comments on Converence Criteria on Page 148. 
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Wing sweep can become a problem for any finite difference scheme employing a 
shearing transformation for the wing planform. The sheared coordinate lines were 
illustrated in Figure 9. As wing sweep increases, the two grid lines collapse on each 
other degrading resolution in the physical Y -direction. Depending on the application, 
the methods performance may be downgraded for wing sweep angles greater than 
50 or 60 degrees. 

Finally, a warning about analyzing complex shapes is in order. Experience 
has shown that the modeling and input of complex geometric shapes is error prone. 

It is hoped that users of this method will familiarize themselves first by analyzing 
simple shapes and slowly progress to more complicated shapes. The array of sample 
cases included in this report should provide a good starting point. In addition, the 
modular construction of the present code (isolated bodies and wings as well as 
wing-bodies) should prove to be useful in diagnosing both code usage and aerodynamic 
problems . 

Future work will include the development of methodology for treating multiple 
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Figure 61. NASA-Douglas Supercritical Wing Transport Configuration 




WING PRESSURE DISTRIBUTION - UPPER SURFACE ~ WING PRESSURE DISTRIBUTION - LOWER SURFACE 

MflCH = 0,750 ALPHA * 2.20 RE = 5.00 M MACH = 0.750 ALPHA = 2,20 RE = 5.00 M 

WING CL =0.381 CM = -0.108 CO = 0.0256 WING CL =0.361 CM = -0.108 CD = 0.0256 


Figure 62. Douglas Transport Configuration Superimposed Computed Pressure Distributions 
= 0.75, a = 2.2° 
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DOUGLAS WING-CYLINDER DOUGLAS WIN6-C IL INDER 

WING STATION 6 2T/B = 0.234 WING STATION II 2T/B = 0.488 

MACH = 0.750 ALPHA = 2.20 MACH = 0.750 ALPHA = 2.20 

5ECTI0N CL = 0.460 CM = -0.061 CO = 0.0340 SECTION CL - 0.554 CM = -0.117 CO = 0.0195 


CD 



DOUGLAS WING-CYLINOER 

WING STATION 17 2T/B = 0.874 

MACH = 0.750 ALPHA = 2.20 

SECTION CL = 0.498 CM = -0.097 CO - 0.0081 

Figure 63. Douglas Transport Configuration Wing Pressure Distribution Correlation, 
IVL =0.75 a = 2.2° 
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DOUGLAS WING-CYLINDER 

WING STATION 11 2T/B = 0.488 

MACH = 0.840 ALPHA =1.85 

SECTION CL = 0.629 CM = -0.145 CD = 0.0210 



DOUGLAS W1NG-CTLINDER 

WING STATION 17 2T/B = 0.874 

MACH = 0.840 ALPHA « 1.85 

SECTION CL * 0.590 CM = -0.105 CD = 0.0012 


Figure 64. Douglas Transport Configuration Wing Pressure Distribution Correlation, 
M 00 = 0.84, 0 = 1.85° 
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DESCRIPTION OF COMPUTER CODE WIBCO 
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GENERAL COMPUTER CODE DESCRIPTION 

The computer code is operational on both IBM and CDC type computers. Over- 
lay structures are not used although this approach (for reducing core requirements) 
may be advantageous depending on facility charging algorithms. The IBM version 
using the extended H compiler (opt = 2) requires approximately 670K 1Q for storage 
and execution. The CDC version with OPT = 1 compiler optimization and segmented 
loader requires about 236Kg for storage and execution. There is considerable use of 
temporary disk storage units. Since interpolation and searching is required, a result 
of the mesh embedding approach, it is useful to have planar potential (v>) arrays 
separate and addressable. As a result, 72 different units (16 CDC units) are cur- 
rently employed. The disk unit number and a description of contents are listed below. 

DISK NUMBER DESCRIPTION 

1 Input data transferred to Unit 1 

7 Flowfield potentials for save/re-start capability 

8 Quick-geometry problem diagnosis printed -output 

10-35 (1 CDC Unit) Global crude grid potential array 

41-58 (1 CDC Unit) Fine wing grid potential array 

61-76 (1 CDC Unit) Fine body grid potential array 

Crude grid wing upper/lower surface boundary conditions 
Fine wing grid upper/lower surface boundary conditions 
Fine wing grid x-coordinate array 
Fine wing grid section surface ordinates 
Crude grid body surface normal (direction cosines) 

Fine grid body surface normal (direction cosines) 

Fine wing grid shearing angles 
Wing and body pressure coefficient arrays 
Boundary layer displacement thickness slope 


80 

81 

82 

83 

84 

85 

86 

87 

88 
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Computer running time will of course vary depending on the facility and the mode 
of operation or operating system. The absolute levels specified may be out dated 
shortly after they are specified. IBM 370 running times are specified below, however, 
since the relative increment for various options will remain essentially steady and 
these increments will be useful in estimating the time and cost of using different 
options . 

TIME (CPU Minutes) 

17 (50 crude/50 crude-fine iterations) 

37 (100 crude/SO crude-fine iterations) 

43 (100 crude/30 crude-fine iterations) 

37 (100 crude/SO crude-fine iterations) 

grid only) 

Wing-body (body modelled in fine grid) 50 (100 crude/80 crude-fine iterations) 

Geometry/Grid verification 1 (No iterations) 

An effort has been made to minimize the amount of data required to define the 
configuration geometry and flow condition. This should simplify matters for most 
applications involving configuration analysis and reduce the chances for input errors. 
For example, the computational grid systems (extent and density) have been set in the 
FORTRAN coding to provide good results under most conditions. Occasionally, it 
will be advantageous to manipulate the preset values and limiters. FORTRAN coding 
changes will be necessary if this is the case. The following values and limiters may 
be modified in certain special applications. 

1) Gas constant ( 7 = 1.4) 

2) Fine wing/body embedded grid limits or extent 


CASE 

Isolated body 
Isolated wing 

Isolated wing w/ viscous interaction 
Wing-body (body modelled in crude 
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3) Fine wing/body embedded grid density 

4) Subsonic relaxation factor ( cj = 1.5) 

5) Boundary layer transition ( X /C tran = 0.05) 

6) The number of inviscid cycles between each viscous calculation 
(currently set to 20) 

All of the sample cases were computed using the basic code without modification. If 
code modifications of type (2) or (3) are made care must be taken to insure that common 
and dimensioned arrays are sufficient and consistent. 

The input data format description can be found on the following pages. The 
descriptions are thought to be relatively straight forward except in the case of wing 
section definition. It is important to extend wing planform/section definition to the 
symmetry plane even for wing -body configurations (see Figure 25). This serves 
several purposes. First, the code will compute a wing -body juncture which will be 
a function of both the configuration geometry and the computational grid system. If 
the computational juncture is slightly inboard of the geometric juncture, section 
definition in this region becomes important. Second, the input planform shape provides 
both the aspect ratio for the lift-induced drag computation and the reference lengths 
and areas used to reduce the integrated pressures to give force and moment co- 
efficients. Finally, the data input for defining a wing -body configuration can be used 
directly for the isolated wing case. This feature can be used to study wing -body 
interference effects. 

It is important to note that there is no provision for input of reference lengths 
and areas. These reference values are computed from the input geometry and printed 
at the end of the output stream. If reference values other than those computed are 
used experimentally, the computed force and moment coefficients must be appropriately 
scaled. 
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LOCATION OF INPUT DATA READ STATEMENTS 


CARD(S) 

ROUTINE 

1-A THRU 3-B 

MAIN 

4-B 

QWIKDE 

5-B THRU 9-B 

CSMDEF 

10-B THRU 13-B 

BLMDEF 


102 



INPUT DATA FORMAT 


Excluding literal cards, all input data cards are punched in seven field ten digit 
format (7F10.0). * A decimal point is required in each field. 


CARD 

NUMBER 

CARD 

COLUMN 

VARIABLE 

NAME 

DESCRIPTION 

Card 1-A 

1-80 

TITLE 

Configuration or run title to identify graphic 
and printed output. 

Card 2-A 

1-10 

CASE 

CASE = 1, Isolated Body (omit cards -W) 
CASE = 2. Isolated Wing (omit cards -B) 
CASE = 3. Wing-Body 


11-20 

AMACH 

Mach Number (AMACH < 1.0) 


21-30 

AOA 

Angle -of- Attack (degrees) 


31-40 

RE 

Reynolds Number (X10^) 


41-50 

AXIT 

Number of initial crude grid iterations. 
(AXIT = 0. for geometry verification only) 


51-60 

AXITF 

Number of crude/fine grid iteration cycles. 


61-70 

VISMOD 

VISMOD = 1. No viscous effects. 

VISMOD = 2. Viscous effects computed at 
end of inviscid analysis. 
VISMOD =3. Inviscid/viscous interaction. 

NOTE: 

: Omit card set 1-W through 5-W for CASE = 1. 

Card 1-W 

1-10 

ASECT 

Number of streamwise sections defining wing 
planform (2 < ASECT < 20). 


11-20 

ANIN 

Number of ordinates defining each wing section 
(ANIN < 60) . 


*NOTE: Card columns 71-80 on card 2-A have been reserved for the variable START 
which provides a solution re-start capability. 


START = 0. 
START = 1. 
START = 2. 
START = 3. 


Conventional solution 
Writes/saves flow field on tape unit 7 

Restarts with tape 7 as input - does not save resulting solution 
Same as 2. But final solution is saved on tape unit 7 
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CARD 

NUMBER 

CARD 

COLUMN 

VARIABLE 

NAME 

DESCRIPTION 

Card 1-W 
(cont'd) 

21-30 

ANOSW 

ANOSW = 0. Sharp nose wing sections. 
ANOSW = 1. Blunt nose wing sections. 


31-40 

XMOM 

X-position about which moment is to be 
computed . 


41-50 

ZWING 

Z -position of wing (waterline). 


51-60 

WPO 

WPO =0. No crude grid output. 

WPO = 1. Crude grid output for diagnostic 
purposes. 


61-70 

ws 

Wing Cp distribution plot scaling per inch 
(typically 0.4 or 0.8). 

NOTE 

: Card set 2 

-W through 5- 

-W is repeated ASECT times. 

Card 2-W 

1-10 

XPL 

Wing section leading edge (X -value) . 


11-20 

YP 

Wing section span position (Y -value). First 
Y -value must be 0,0 (symmetry plane) , even 
for wing -body case. 


21-30 

XPT 

Wing section trailing edge (X-value) . 


31-40 

TWIST 

Wing section local incidence (twist angle in 
degrees) . 


41-50 

AKODE 

AKODE = 0. Section ordinates identical to 
preceding section (omit cards 
3-W through 5-W). 

AKODE = 1. New section definition expected 
on cards 4W and 5W 

Card 3-W 

1-70 

XINW 

Wing section X -coordinates (cards 3-W 
defined only for first wing section, ANIN 
values expected) . 

Card 4-W 

1-70 

YINU 

Wing section upper surface Y-coordinates 
(ANIN values). 

Card 5-W 

1-70 

YINL 

Wing section lower surface Y-coordinates 
(ANIN values). 


NOTE: Omit card set 1-B through 13-B for CASE - 2 


Card 1-B 


1-10 


BKOD = 1. 
BKOD = -1. 

BKOD = 2 . 

BKOD = -2 . 


Infinite cylinder (only RADIUS need be input) . 

Same as BKOD - 1. No embedded body grid. 
Crude grid body representation only . 

Simple axisymmetric body definition requested 
(input BXIN, RIN on card(s) 2-B and 3-B). 

Same as BKOD = 2 . No embedded body grid . 
Crude grid body representation only. 
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CARD 

CARD 

variable 



NUMBER 

COLUMN 

NAME 

DESCRIPTION 

Card 1-B 
(cont T d) 


BKOD = 3 . 
BKOD = -3. 

Complex body definition requested (input Quick- 
Geometry model on card(s) 4-B through 13 -B). 

Same as BKOD = 3. No embedded body grid. 
Crude grid body representation only. 


11-20 

21-30 

BNOSE 

BTAIL 

Body nose (X-value)| 
Body tail (X -value) j 

For BKOD = ±2. or ±3. 


31-40 

BNIN 

Number of axisymmetric body coordinates to 
be input. BNIN < 60 (for BKOD = ±2 only). 


41-50 

RADIUS 

Cylinder radius for BKOD = ±1 only. 


51-60 

ANOSB 

ANOSB = 0. Sharp nose body ) n 
ANOSB = 1. Blunt nose body J BKOD = ±2 onl y • 


61-70 

BS 

Body Cp plot scaling per inch (typically 0.08). 

NOTE: 

Omit card 

sets 2-B and 3 

-B for BKOD = 1 or BKOD = 3 . 

Card(s) 2-B 

1-70 

XINB 

Axisymmetric body X-coordinates (BNIN 
values). 

Card(s) 3-B 

1-70 

RIN 

Axisymmetric body radii (BNIN values) . 

NOTE : 

Omit card sets 4-B through 13 -B for BKOD = 1 or BKOD = 2. 

Card 4-B 

1-70 

VTITLE 

Quick-Geometry model title. 

Card 5-B 

1-10 

ACSM 

Number of distinct cross-section models 
(ACSM card sets 6-B and 7-B will follow) . 

Card 6-B 

1-10 

ADUM 

Running count of current cross-section model 
(1-ACSM) . 


11-20 

AARC 

Number of arcs in current cross-section 
model (AARC Card(s) 7-B will follow). 


21-60 

CTITLE 

Title or descriptor of current cross-section 
model. 

Card 7-B 

1-8 

arcnam 

Arc or component name. 


11-14 

ASHAPE 

Arc or component shape. 


21-28 

PNTNAM(l) 

Control point name for beginning of this arc. 


31-38 

PNTNAM(2) 

Control point name for termination of this arc . 


41-48 

PNTNAM(3) 

Slope control point name for this arc, if 
required. 

Card 8-B 

1 

h-* 

O 

ANTCSM 

Number of cross-section models to define 
entire body (ANTCSM card(s) 9-B will 
follow) . 


105 


CARD 

NUMBER 

CARD 

COLUMN 

VARIABLE 

NAME 

DESCRIPTION 

Card 9-B 

1-10 

ADUM 

Running count of current cross-section model 
(1-ANTCSM) 


11-20 

AMODEL 

Index corresponding to already defined cross- 
section models (between 1 and ACSM). 


21-30 

XCSMS1 

Starting X-station for current cross-section 
model. 


31-40 

XCSMS2 

Ending X-station for current cross-section 
model . 

Card 10-B 

1-10 

BLINE 

Number of body line models to be defined by 
segments (BLINE card sets 11-B and 12-B 
follow) . 


11-20 

ALIAS 

Number of body line models to be aliased 
(Alias card(s) 13-B follow). 

Note: 

Card set 11- 

B and 12-B 

is repeated BLINE times. 

Card 11-B 

1-10 

BLSEG 

Number of segment(s) defining body line 
model. 


11 

BYORZ 

The letter Y or Z indicates which data 
definition is to follow. 


12-19 

BNAME 

Body line name to be defined. 

Card 12-B 

1-4 

SSHAPE 

Segment shape. 


11-20 

D(l) 

X-station for beginning of segment. 


21-30 

D(2) 

Y or Z value corresponding to D(l). 


31-40 

D(3) 

X-station for termination of segment. 


41-50 

D(4) 

Y or Z value corresponding to D(3). 


51-60 

D(5) 

X-station for segment slope control point. 


61-70 

D(6) 

Y or Z value corresponding to D(5) 

Note: 

Card set 13 

-B is repeated ALIAS times. 

Card 13-B 

11 

BYORZ 

The letter Y or Z indicates which data 
definition is to follow. 


12-19 

BNAME 

Body line name to be defined. 


21 

AYORZ 

The letter Y or Z indicates which definition 
is to be used for aliasing. 


22-29 

ANAME 

Body line name to which BNAME is aliased.) 
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SAMPLE INPUT DATA SETS 
(8 CASES) 


NOTE: These data sets were used for the correlation studies 
(See Comparisons Section, page 55) 
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Case No. 1 

NACA RM L53F07 ISOLATED BODY 


1 . 0 

O.on 

8. 4 

1.0 

50. 

5 0. 

1.0 

7 . 

0. 

. . 1023 

32. 

0 . 

0 . 

on 

O 

« 

n i 

0. 

0.5 

0. 75 

1.25 

7.50 

5.00 

7.50 

10.0 

1 5.0 

70.0 

25. 0 

30.0 

35.0 

4 0.0 

4 5.0 

50.0 

55.0 

60. 0 

65. 0 

70.0 

75.0 

00.0 


89.0 

90.0 

91.0 

92. 0 

93.0 

94.0 

O 6. 0 

96. 0 

ino.o 




0 . 

0,731 

0, 798 

O. 428 

0.722 

1.205 

1. 613 

1.071 

2. 59 3 

3. 090 

3.465 

3. 74 1 

3.933 

4.063 

4.143 

4.167 

4. 130 

4.024 

3.842 

3.562 

3. 128 

2,526 

1.952 

1. 265 

1.1 20 

1.005 

0.925 

0.872 

0.041 

0 . 03 9 

0. 833 

0.833 








Case No. 2 






NACA TND 

7331 ISOLATED BODY 



1 .0 

0.99 1 

0 . 0 

1.0 

50. 

5 0. 

1.0 

2. 

0 . 

100. 

44. 

0 . 

1. 

3.03 

0 . 

0.2 

0. 4 

0.6 

0.8 

1.0 

1 .2 

1.4 

1 .6 

1.8 

2.0 

4.0 

6.0 

3.0 

10.0 

12.0 

14.0 

16. 0 

18.0 

20. 0 

24.0 

2 8.0 

32.0 

40.0 

48. 0 

50. 0 

52. 0 

50.0 

68.0 

7 2.0 

76.0 

80. 0 

82.0 

84,0 

86.0 

98.0 

90.0 

92 .0 

94. 0 

96.0 

97. 0 

98.0 

99.0 

100. 0 






0 . 

0.696 

0. 946 

1.134 

1.291 

1.422 

1. 53 9 

1 .645 

1.742 

1. 832 

1.916 

2.556 

3.005 

3.357 

3.648 

3.994 

4. 108 

4.295 

4. 460 

4.607 

4. 855 

5.052 

5 . 207 

8. 412 

c .497 

5.500 

5.497 

5. 412 

5.207 

5.052 

4. 855 

4.607 

4. 460 

4.295 

4 . 103 

3.894 

3.649 

3. 357 

3.005 

2.556 

2. 271 

1.916 


1.422 1.427 
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Case No. 3 



GULFSTREAM II 

ISOLATED 

FUSELAGE 

NOSE 

1 • 

o. ao 

3.1 

1*0 

5 0. 

50. 

3 » 

A* 0 

200. 0 

0. 

0. 

0. 

GULFSTREAM FUSELAGE QUI CK- GEOM ETR Y MODEL 


3* 






1 • 

2* 

NOSE TO 

W INDSH IELD 

BASE 


BCDYLO 

ELLI 

8DYBCL 

BDY MHB 

BOYL SCP 


BODYHI 

ELLI 

0DYMHB 

BOYTCL 

BDYUSCP 


2. 

3* 

W INDSHIELD 



CANLO 

ELLI 

0DYBCL 

BDY M HB 

bdylscp 


CANS I 

ELLI 

9DYMHB 

CANLOW 

C ANLSCP 


WINDF 

LINE 

C ANLOW 

BDY T CL 



3* 

4* 

CANOPY 




C ANCPLC 

ELL I 

BDY0CL 

BDY MHB 

BDYL SCP 


CANOPS I 

ELLI 

BDYM hB 

CANLOW 

C ANLSCP 


W INDSI 

LI NE 

C ANL 0 ta 

CANH IE 



CANOPUP 

4* 

ELLI 

C ANN I E 

BO Y T CL 

CANT SCP 


1 • 

1 • 

8.0 

64. 0 



2. 

2* 

64.0 

83. 5 



3. 

3. 

83.5 

l 33. 0 



4* 

1* 

133.0 

200.0 



10* 

1 1 • 





2* 

ZBDYBCL 





ELL X 

8.0 

78.5 

133. 0 

53.0 

8.0 

LINE 

1 33.0 

53.0 

200. 0 

53.0 


4* 

ZBOYTCL 





ELLX 

8* 0 

78.5 

64.0 

113.5 

a. o 

LI NE 

64* 0 

1 13.5 

83.5 

1 32. 5 


CUB I 

83*5 

1 32. 5 

133.0 

1 47.0 

9 8. 0 

L INE 

133.0 

147. 0 

200.0 

147.0 

1 • 

YCENTER 





LINE 

8.0 

0 .0 

200.0 

0.0 


2. 

ZBDYMHB 





CUB! 

8. 0 

78.5 

18 1.0 

1 00. 0 

1 35.5 

LINE 

181.0 

100.0 

200.0 

100.0 


1 • 

YBDYMHB 





ELL X 

8.0 

0. 

200. 0 

4 7.0 

8.0 

1 • 

ZC ANLOW 





L f NE 

64.0 

113.5 

133. 0 

113.5 


1* 

ZC ANHI E 





L I NE 

83.5 

1 32.5 

133.0 

132.5 


2* 

YC ANLOW 





L INE 

64.0 

0. 

79.5 

26.0 


L INE 

79.5 

28.0 

133. 0 

42.0 


2* 

YC AN HIE 





LI NE 

83.5 

0 . 

96.0 

24. 5 


LI NE 

96.0 

24.5 

133.0 

34. 0 


1* 

YC ANT SC P 





L INE 

83.5 

0. 

133.0 

22.0 



ZMAPAX I S 

Z 8DY MHB 





YMAPAX I S 

YCENTER 





ZCANTSCP 

ZBOYTCL 





Z C ANL SCP 

ZC ANLOW 





YC ANLSCP 

YBDYMHB 





YBDYTCL 

YCENTER 





YBOYBCL 

YCENTER 





YBCYL SCP 

YBDYMHB 





Z8DYL SCP 

ZBDYBCL 





Y8DYUSCP 

YBDYMHB 





ZBOYUSCP 

Z0DY T CL 





i *o 

0*4 


53* 0 

94*5 

147.0 

1 00*0 
4 7* 0 


sd ^ O \>J tvj 


Case No. 4 

ONERA M6 ISOLATED WING 


♦ 

*56354 
47. «1 197 
90.6199S 

0 . 

3.57742 
4.7166 1 
1.22399 

0 . 

-3. 57742 
-4.7166 1 
- 1.2238° 
0.85670 


2 . 

2 _* 

0. 

0 . 

10. 859 
47.644 
9 1.701 

0. 

3.946 

5.916 

2.297 

3. 

-2.939 

-3.969 

-1.903 

45.625 


0 . 9U 07 

3.01 

2.5 

100. 

• 

O 

CD 

1 . 0 

2 4. 

. _ 

0.452 

0. 

1, 

0.4 

0. 

1.0 

0. 

0. 



0. 29441 

0. 457 04 

0.084 13 

1. 94609 

2.97008 

5. 41248 

1 6. 4 q Q76 

21 .870 q 6 

27. 17978 

32. 42726 

37. 61446 

42.74223 

9 6. 30° 3 7 

60.23757 

65. 11093 

69.0 3027 

77. 05998 

84.07324 

96. 6 1 96 0 

100.0 





n # 97 q 59 

1. 17419 

1.68984 

2. 24 c 45 

2.58245 

3 . 03278 

4.1 9099 

4. 50507 

4.71QR7 

4.84 902 

4.89296 

4 . 84 8 88 

4. 36 74 1 

4. 0524 1 

3.68990 

3. 29 402 

2.65505 

1.95838 

0.48907 

0. 07052 





-0. «7069 

- 1. 17 4 19 

- 1. 68984 

- 2. 24545 

- 2. 58245 

-3.03278 

-4. 1 9039 

-4.50507 

-4. 71 qp7 

-4. 8 490 2 

-4. 99296 

-4. 84988 

- 4. 3 6 74 1 

-4. 05241 

-3. 6 9 90 0 

- 3. 2 5 402 

- 2. 65505 

- 1.95838 

-0. 4 890 7 

-0.07052 





1. 49 39 

1. 3944 

c. 

0. 






Case No. 5 





NACA RM 

A9K01 ISOLATED WING 



0 .8500 

4. 00 

2.0 

100. 

80. 

3 . 

26, _ ... 

!*.0 

25. .4 80 

0. 

1,0 

0.6 

0. 

26 .67 

0. 

0. 


8. 129 

0. 465 

0. 7 33 

1.275 

2. 644 

5.388 

1 6.279 

21 . 64 7 

26. 959 

32.213 

37. 413 

42.555 

6 2. 6 0 4 

57.649 

62. 569 

67.4 33 

72.242 

76.998 

86. 3 50 

90.948 

95. 497 

100. 0 



0.00 8 

1. 103 

1.411 

1. 96 1 

2.754 

3.355 

4.614 

5. 175 

5.580 

5. 84 5 

5.978 

5. 983 

5.525 

5. 135 

4.666 

4. 133 

3.55 1 

2.934 

1.662 

1. 049 

0.4R4 

0. 



-0. 8 20 

-0.970 

- 1. 221 

- 1.632 

- 2. 1 96 

-2.608 

-3.4 90 

-3.704 

-4. 03 5 

-4. 177 

-4. 220 

-4.155 

- 3. 67 3 

-3. 307 

-2. 887 

-2.432 

- 1. 954 

- 1.471 

-0.573 

-0.21 6 

0.022 

0. 



6 0. 0 

68. 955 

-0.6 

0. 
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3 • 

2* 

12*6 25 

3 • 

10*0 

4 5*3 
30*0 

3 • 

1 • 824 
2.99 2 
1.43/ 

3. 

- 4 . 8 24 
-2*992 

- I .4 37 

2 5.3 75 

2 . 

3* 

4* 0 
18*0 
32*3 

3* 

0* 78 84 
1 *65 72 
1*0104 


Case No. 6 

NACA RM L51F0 7 AXI-SYM BODY/45 DEG SWEPT WING 


80 . 2*0 

uo 0.6 


0.93 

4.0 

26 o 

1.0 

0 * 

20.125 

3.5 

0.75 

1 5 o 9 

20.0 

50.0 

55*0 

85.0 

9 0.0 

0 * 464 

9 * 56 3 

2.194 

2 . 474 

2.925 

2. 793 

l . 033 

0 * 727 

- 0.464 

- 0*563 

- 2.194 

- 2.474 

- 2 * 925 

- 2. 793 

- 1 .033 

- 0.727 

12.0 

29 . 8'5 

0 a 

33.333 

0.2 

0.3 

6.0 

6.0 

20 . 0 

22 . 0 

33.0 

33 . 333 

3.0924 

3.1192 

1 .0372 

1 o 2360 

1 .6668 

1 .6520 

Oo 88 

0.833 


2 o 0 2 

l 00 . 

20 . 

0.0 

3 . 

1 . 0 

1 .25 

2*5 

25 . 0 

3 9 . 0 

60.0 

65 * 0 

95.0 

1 00 . 0 

0.718 

3*93 1 

2.637 

2 . 842 

2 . 6 0 2 

2.364 

0.370 

3 . 01 3 

- 0 . 71 8 

- 0*981 

- 2.687 

- 2.342 

- 2. 602 

- 2*364 

- 0.37 0 

- 0*013 

0 . 

3 . 

25 . 

0 . 

0.5 

1.0 

10.0 

l 2.0 

24.0 

26.0 

33.8 


0.1712 

0.2688 

1 • 3860 

1 . 4964 

1.6096 

1*5368 

0.833 



5 . 0 

7.5 

35.0 

40 . 0 

70 . 0 

75.0 

1.313 

1.591 

2 . 945 

2. 996 

2 . 087 

1 a 775 

- 1.313 

- 1*591 

- 2.945 

- 2.996 

- 2 . 087 

- 1. 775 


0*0 

0.4 

2.0 

3.0 

14 . 0 

16.0 

28 * 0 

30 . 0 

0.4823 

0 . 6452 

1*5732 

1 • 6252 

1*4248 

1.2512 
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3 9 

6 • 

2 Q .66 

0 . 

10 * 

45 . 

80 * 
3*740 
6*572 
8 * 76 6 
3 * 44 5 
0*00 0 
- 3 * 6 6 b 
- 6*0 44 
- 2. 3 09 
32*0 2 
30*38 
0 * 53 6 
5 * 51 l 
7.72 5 
3*188 
0.000 
- 2 * 6 46 
- 5.0 07 

- 1 .998 
41*4 1 
0.37 4 
4 • 69 9 
6.736 
2. 929 
0.000 

- 1 .7 87 
- 3.782 

- I .5 93 
50.8 1 
0. 30 6 
4. 42 3 
6 . 370 
2* 90 4 
0.00 0 
-1.187 

- 2 • 5 63 

- 1 . 0 39 
69.2 7 
- 2 . 

0 . 

1 4,0 
70.0 
108. 3 

0 . 

4*95 

5. 04 
2. 60 


BOEING 


• 78 

2.0 

26. 

0* 

0.0 

70.27 

0*5 

0* 7 5 

15. 

20. 

50. 

55 • 

8 5# 

90. 

1 • 750 

2.055 

7* 66 i 

3.36 0 

3.258 

7. 62 0 

2.564 

1 • 723 

-0.872 

-1.076 

-4.51 0 

- 5. 228 

-5.678 

-5. 185 

- 1 . 732 

-1.155 

4*3 7 

62* 89 

1 0. 36 

67. 06 

1 *399 

1*654 

6* 436 

7* 07 0 

7m 362 

6. 895 

2*391 

1. 594 

-0. 733 

-0. 895 

-3*318 

-3.925 

-4.793 

- 4 * 4 48 

- 1 .498 

-0*999 

l 6. 73 

59*20 

l . 123 

1.3^0 

5 . 530 

6. 047 

6. 5 09 

b* 1 36 

2*197 

1*465 

-0.6/0 

-0 .684 

-2* 265 

-2.731 

-2.696 

-3*487 

-1.194 

- 0.796 

29* 30 

6 4.90 

0*997 

1.198 

5 . 2 l 0 

5. 704 

6.178 

5. 850 

2.178 

1.452 

-0.465 

- 0 . 558 

- 1. 4 79 

-1.770 

-2*50 1 

-2. 330 

-0.775 

-0.520 

53* 96 

74* 97 

0.C 0 

108.30 

2.0 

4.0 

16* 0 

20. 0 

78.0 

80 • 0 

2*50 

3. 35 

5*04 

5. 04 

5*04 

4.92 


Case No. 7 
KC-135 TRANSPORT 


5 . 25 

1 00 * 

52 . 90 

- 4 . 25 

1.8 

l . 0 

1*25 

2.50 

25 . 

30 . 

60 . 

65 . 

95 . 

1 00 . 

2.539 

3*47 0 

8.860 

9.137 

6 . 870 

6 * 029 

0.86 1 

0 * 

- 1 .390 

- 1.9 15 

- 5 . 757 

- 6. 1 04 

- 4.618 

1 

. 

O 

>-• 

- 0 . 577 

0 . 

1.0 

0 * 

1.76 

1 • 

2.060 

2*840 

7.500 

7 * 76 0 

6.284 

5*570 

0 . 797 

0 * 

- 1.123 

- 1 • 4 75 

O 

* 

1 

- 4 * 746 

- 3 . 995 

- 3*496 

- 0 * 499 

0 . 

1.7 0 

1 *0 

1.687 

2*350 

6*422 

6*666 

5 . 643 

5 * 049 

0 . 732 

0* 

- 0. 837 

- 1 * 053 

- 3 . 1 1 8 

- 3 * 425 

- 3 * 174 

- 2 * 787 

- 0. 398 

0 . 

1 .44 

1 • 0 

1.519 

2*160 

6 . 050 

6.280 

5 . 41 2 

4 * 880 

0 . 726 

0 * 

- 0.657 

- 0 * 750 

- 2.059 

- 2*275 

- 2.078 

- 1*018 

- 0.260 

0 . 

0 . 80 

0 * 

22 . 

5.04 

6.0 

8.0 

30 . 0 

40 . 0 

o 

• 

10 

CO 

90 . 0 

3.90 

4 . 30 

5.04 

5.04 

4.65 

4.20 


00 * 

1 • 

1 * 

0.8 

5.0 

7.5 

35 . 

40 . 

70 . 

75 . 

4.830 

5.813 

9.218 

9.105 

5 . 168 

4*300 

- 2 . 61 3 

- 3 . 178 

- 6.275 

- 6 * 250 

- 3. 464 

- 2*886 


3.998 

4 .850 

7*920 

7*9 10 

4.782 

3*985 

- 1*915 

- 2*290 

- 4. 980 

- 5*070 

- 2*996 

- 2.497 


3*360 

4.114 

6 * 805 

6.829 

4 * 380 

3*687 

- 1*307 

- 1 * 544 

- 3*642 

- 3 . 763 

- 2*389 

- 1.991 


3*130 

3.857 

6.430 

6.449 

4.267 

3.605 

- 0.895 

- 1.041 

- 2 * 445 

- 2.543 

- 1*558 

- 1*299 


1 . 

• 4 

10.0 

12.0 

50 . 0 

60.0 

95.0 

100*0 

4*60 

4*30 

5.04 

5*04 

3*55 

2 .65 
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Case No. 8 

NASA TM X-3431 TRANSPORT CONFIG. 


3. 

0.909 

3. 91 

1 8. 

2 5. 

1. 0 

36. 3 

0 . 

73.296 

0 . 

0.00 25 

0. C05 0 

0.0750 

0.100 

0. 150 

0.500 

0.550 

0. 600 

0.850 

0.900 

0. 950 

-0.0225 

-0.0 153 

-0.0119 

0.0140 

0.0182 

0. 0242 

0.0173 

0.0124 

0. C 07 4 

-0.0212 

-0.0279 

-0.0343 

-0.0225 

-0. 0 296 

-0.0321 

-0.0618 

-0. 0664 

-0.0728 

-0.0797 

-0.0751 

-0.0700 

-0.0434 

-0.0424 

-0.0435 

4 5.5 

5.96 1 

'8.407 

53.2 

1 1.041 

79 . 27 ^ 

-0.0502 

-0. 0 417 

-0.0368 

-0.0110 

-0. 0074 

-0.00 4 0 

-0.0083 

-0. 0 110 

-0.01 44 

-0.0384 

-0.0440 

-0.0505 

-0.0502 

-0. 0 574 

-0.0607 

-0.0910 

-0.094 9 

-0.10 00 

-0.0988 

-0.0954 

-0.0908 

-0. 0631 

-0. Or, 1 1 

-0.06 11 

58.662 

14. 7 23 

60. 426 

-0.0880 

- 0.079R 

-0.0^52 

-0.05 1 1 

-0.0 479 

-0.0434 

-0.0398 

-0.0420 

-0 .0446 

-0.0624 

-0. 0672 

-0.0735 

-0.0880 

-0. 9 959 

-0 .0995 

-0. 1263 

-0. 1297 

-0. 1332 

-0. 1296 

-0. 1 261 

-0.1214 

-0.0882 

-0.0857 

-0.0862 

6 3. 248 

1 7. 779 

81 .748 

-0. 1135 

-0. 1056 

-0. 1021 

-0.0772 

- 0 . 0 73 4 

-0.06 95 

-0.0625 

-0. 0 638 

-0.0652 

-0.0791 

-0.083 6 

-0,0900 

-0. 1 1.15 

-0. 1 20 9 

-0. 1243 

-0.1482 

-0. 1 51 4 

-0.1552 

-0. 1504 

-0. 1 467 
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OUTPUT DATA (PRINT AND PLOT) FORMAT 

Printed and plotted output data is provided. Since a typical print or plot sequence 
is lengthy, only a brief description of each type of output will be provided here. Note 
that samples of the plotted output can be found in the results and geometry verification 
sections of this report. 

Printed Output 

The printed output can be divided into three distinct sections. 

Section I Input Data Listing/Geometry and Grid System Verification 
Section II Relaxation Solution Convergence History 

Section III Computed Velocities, Pressures, Forces, Moments, Reference 
Lengths and Areas 

Within each section, the output data will be printed in the following sequence. 

Section I 

• Input Data Listing 

• Case Flow Condition 

• Nominal Extent of Fine Embedded Grid Systems 

• Quick-Geometry Model Error Diagnostic Information (BKOD = ±3 only) 

• Configuration Position in Global Crude Grid System 

• Body Crude/Fine Grid Limiters 

• Global Crude Grid Coordinates 

• Fine Embedded Wing Grid Coordinates 

• Fine Embedded Body Grid Coordinates 

Section II 

Phase 1 (Global Crude Grid Solution) 

The following information is printed for each crude grid iteration. 


A 


• Iteration Count 

• Maximum Correction to the Flow Field Potential (A#> 

• Grid Position of A** 

MAX 

• Wing Spanwise Circulation (T) Distribution 


MAX 


) 
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The following information is printed at the end of the Global Crude Grid Solution. 

• Section Mach Chart (0 indicates subsonic flow) 

(1 indicates supersonic flow) 

B J (8 indicates wing section surface) 

(7 indicates wake surface) 

• Span row number/ *1 -position/local chord length, chord position, pressure 
■ coefficient, disturbance velocity, section circulation (T). 

Phase 2 (Crude/Fine Grid Interactions) 

The following information is printed for each crude/fine grid interation cycle. 

• Phase 1 (A) Output for Fine Embedded Wing Grid System 

• Phase 1 (A) Output for Fine Embedded Body Grid System (for BKOD >0 only) 

• Phase 1 (A) Output for Global Crude Grid System 

• Wing Spanwise Circulation (D Distribution 

For viscous interaction cases, the following information is printed every 20th 

cycle . 

• Wing Upper/Lower Surface Boundary Layer Separation Point (x/c) 

• Boundary Layer 5* Slope Added to Wing Boundary Conditions. 

Section III 

The following information is printed at the end of the solution process. 

• Phase 1 (B) Output for Wing Embedded Fine Grid System 

• For (VISMOD =2,3) Wing Upper/Lower Surface Boundary Layer Separation 
Point (x/c) 

• Wing Section C . , C m , C d , C f 

cu OJ OJ OJ 

Note* C is section moment about local quarter chord position, 
m 

OJ 

• Spanwise Load, Moment, Drag Distributions 

• Wing Exposed Area, Total Area, Aspect Ratio, Taper Ratio, Mean Aero- 
dynamic Chord, Average Chord, X-position about which Moments are computed 

• Total Wing C , C , C 

wING 1V WING WING 
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• Wing Pressure Drag C_ 

P(w) 

• Wing Friction Drag C v 

* WING 

• Body Grid Mach Chart 

• Body Angular Cut Pressure/Velocity Distribution 

• Body Longitudinal Load and Drag Distribution 

• Body Length, Wetted Area, Projected Area, Max. Cross-sectional Area, 
Reynolds Number based on body length 

• Total Body C T , C C 

bODY m BODY U P(b) 

• Body Pressure Drag C n 

P(b) 

• Body Friction Drag C^ 

1 B 

• Total Configuration C , C C 

B M D 

• Wing-Body Wave Drag C n 

d wave 

• Wing-Body Friction Drag C 

F 

• Wing-Body Lift-Induced Drag C 

Dl 

• Wing-Body Spanload Efficiency Factor (E) 

Plotted Output 

The plot output can be divided into two separate sections. 

Section I (Input Geometry Verification) 

• Title/Case/Flow Condition Label 

• Body Cross-Sections 

• Input Wing Sections 

• Configuration Planview 

• Configuration Head-On View 

Section II (Computed Results) 

• Superimposed wing pressure distributions (upper/lower surface) with total 
wing C L , C M , C D label. 

• Wing planform with section shapes at computed span stations. 
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Detailed wing section pressure distributions with section C^, C m> C d label. 
Detailed body angular cut pressure distributions. 


Body longitudinal load plot with body C, , C , 

B B 

C D label. 


C and wing -body C T , 

Dt, 


C 


M’ 


Wing -body spanload plot with span efficiency, lift-induced drag, wave drag, 
friction drag label. 



INPUT GEOMETRY VERIFICATION 

Coding for the graphical inspection of input wing and body geometry has been 
included. Since the program requires considerable computer time and core storage 
to operate and some facility budgets may not provide for a number of error -filled 
submittals, it is recommended that the geometry verification mode of operation be 
used before submitting for a complete and expensive relaxation solution. The graphic 
output coupled with printed output for geometry and grid systems should be sufficient 
to diagnose user input errors. In particular, it has been found advantageous to make 
the geometry check and perform a single crude and a single fine grid iteration (very 
cheap) to check code flow before a complete analysis is performed. 

The body cross-section array is first to be plotted after the case title and flow 
conditions are listed (see Figure A-l). In this case, an error in specifying a Z- 
coordinate of the canopy top centerline is apparent. In Figure A-2, an error in 
specifying a Y-coordinate of the windshield base is illustrated. Finally, Figure A-3 
illustrates the corrected and final shape of the fuselage to be analyzed. 

The defining wing chord sections are displayed after the body geometry (see 
Figure A-4). Each is blown up to a ten inch chord so errors in coordinates will 
become visible. It is important to note that there is no mapping involved in the 
present method, thus, no coordinate smoothing or manipulation is employed. As a 
result, what you input is what you get. Irregularities in input coordinates will cause 
oscillations in computed pressure distributions. 

A configuration plan -view will follow the wing section figures. This should be 
used to insure that the wing planform is properly defined and its placement on the 
fuselage is correct (see Figure A -5). Finally, a head-on view is plotted (see Figure 
A-6). Once again, check to see that the wing and body are indeed attached. 
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Figure A-1. Sample Input Geometry Verification Plot with Error in Z-Coordmate 
of Canopy Definition 



Figure A-2. Sample Input Geometry Verification Plot with Error in Y-Coordinate 
of Canopy Definition 
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GULFSTRERM II ISOLATED FUSELAGE 






Figure A-4. Sample Input Geometry Verification Plot for Wing Section! 



PLAN VIEW 


Figure A-5. Sample Input Geometry Verification Plot, Plan View 



Figure A -6. Sample Input Geometry Verification Plot, Head-On View 



SUBROUTINE CALL SEQUENCE 


MAIN 



NOTE: 



Fuselage Geometry Definition 
Fuselage Geometry Interrogation 
Relaxation/Boundary Layer Routines 


SPLINE 

© 

— © 

SAREA 


125 



QWIKDE 


SUBROUTINE CALL SEQUENCE 



CSMDEF 


DSETUP 


DLOKUP 


CSMCHK 


BLMDEF 



BLMCHK 


I GEMOUT 
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SUBROUTINE CALL SEQUENCE ( B 


QWIKLO 


BLMSET 

BLGEOM 


CSGEOM 


CSMSET 


C BLGEOM 
CSMCOE 
I — AORDER 


— THELIM 


CSMFLT 


I — CSCALC 


1 — CSMINT 
— SINCOS 

VDOTV 

— MDOTV 


CSCALC 

BLGEOM 

CSMCOE 

LIN LIN - 

LINELL - 

' ELLELL- 


SINCOS 

VDOTV 

MDOTV 

CROSS 

ELLCAL 

ELLCAL 

CROSS 

ESTNXT 

SETNXT 
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CONTRL 


SUBROUTINE CALL SEQUENCE 



i TRID 


GLOBAL — 


BODYC 


TRID 


bodval 


INTURP 

INTERB 

■POCRUD 

•WINGF 1 


FILL 


TRID 

INTURP 


— BODYF 


— BODFIX 


-POFINE- 


TRID 

-PLOTER 


•INTEG 


SPLINE 


-LED RAG 


MAST 


~® 


■FIT2 

•DELTA1 

-INTRP 


■ SERIE S — RFAST — FAST 


•SPLINE 


■BODFM ■ 


' — BLAR- 


- VINTER 


SPLINE 


L@ 


® 


SMOTH 


NOTE 


: @ 


BOUNDARY LAYER ROUTINES 
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SUBROUTINE CALL SEQUENCE ( D 


STRIPK 


BLLAM 


SLOPY 

SPLN1 

SPLN1X 

I ORDIN 


VBRAD 


VNUSUB 

SLOPBL 

TANCAL 

SPRINT 

FINT 

SOLVEB 

ORDIN 

REDUCX 

RLORD 

GORD 


GRAD 


' SIMPSN 



SUBROUTINE DESCRIPTION 


AORDER 

BLAR 


BLGEOM 

BLLAM 

BLMCHK 

BLMDEF 

BLMSET 

BODFIX 

BODFM 

BODLIM 

BODVAL 

BODYC 

BODYF 

CONTRL 


Orders a set of numbers by permutation index. 

Main control routine for laminar and turbulent modified chordwise 
boundary layer calculation. Computes boundary layer displacement 
thickness (6*) slope for viscous/inviscid interaction mode of operation. 

Assigns body line model values and derivatives to control point 
coordinates. 

Computes Thwaites laminar boundary layer with Rott and Crabtree 
compressibility modification. 

Correlates and checks the input data deck and the indices for the 
generated body line math models. 

Defines body line models from the input data. 

Controls the determination values and first and second derivatives for 
all body line models at a given x-station. 

Computes potentials on fixed wing/wake surface in body fine grid given 
solution in fine wing grid and global crude grid. 

Computes integrated body force and moment coefficients. 

Computes J and K limiters for body boundary in both crude and fine 
grid systems. 

Computes body boundary point potential values. 

Finite difference approximations and relaxation solution for body 
boundary in crude grid. 

Finite difference approximations and relaxation solution for fine body 
grid. 

Main control routine for relaxation solution of governing equation, 
interpolation, boundary layer analysis and printed/plotted output. 
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CROSS 

CSCALC 

CSGEOM 


CSMCHK 

CSMCOE 

CSMDEF 

CSMFLT 

CSMINT 

CSMSET 

CURVES 

DELTA1 

DLOKUP 

DSETUP 


ELLCAL 

ELLELL 


Solves for the interaction of two lines in a plane. 

Computes radial position and derivatives for specified cross-section 
model, arc, and 0 

Is the main subroutine in the look-up portion of the QUICK system. 

It is called to establish r' = f (0’, x). It calls appropriate subroutines 
to evaluate body line values and construct cross-section geometry at 
a given x-station. It is used for all geometry model interrogation. 

Correlates and checks the input data deck and the indices for the cross- 
sectional math model. 

Composes the equations which are to define the cross-section geometry 
at a given station. 

Logically defines the cross-section models from the input data. 

Creates control point definitions to permit the insertion of a smooth 
fillet between cross-sectional arcs. 

Locates user specified intersections between cross-sectional arcs 
and adjusts their use-theta limits. 

Sets up the control point coordinate arrays used to define the cross- 
section geometry at a specified x-station. 

Calculates values and first and second derivatives for individual curve 
fits. 

Interpolation routine for wing spanload. 

Is a simple dictionary look-up routine. It assigns an index to match 
an input name to a codeword list, but is not capable of adding new 
items to that list. 

Is an adapting dictionary look-up routine. New items are added to a 
codeword list, an index (counter) is returned for the codeword, and an 
indicator (INEW) is set equal to 1 when a new item is encountered. 

Set up for ellipse. 

Calculates intersection of two ellipses. 
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ESTNXT 

FAST 

FILL 

FINT 

FIT2 

GEMOUT 

GEVER 

GLOBAL 

GORD 

GRAD 

INTEG 

INTERB 

INTRP 

INTURP 


KRVDEF 

LID RAG 

LINELL 

LINLIN 

MAIN 

MAST 


Estimates non-linear root by modified inverse quadratic. 

Fast Fourier transform of complex data. 

Performs interpolation controlled by INTURP. 

Simultaneous triple interpolation. 

Determines cubic spline fit coefficients for input spanload distribution. 

Ensures that all body lines required by a cross-sectional model are 
defined for the range of that model. 

Controls geometry verification plotting. 

Finite difference approximations and relaxation solution for global 
crude grid. 

Bradshaw's G function. 

Slope of a function at its tabulated points . 

Integrates wing load distributions for lift, moment and drag 
coefficients . 

Interpolation routine for body fine/global crude grid communication. 

Converts input spanload distribution to a fine over spaced distribution. 

Controls interpolation for filling fine mesh points using crude grid 
potential values. Updates crude mesh given fine solution. Updates 
fine mesh given crude solution. 

Calculates coefficients for the various curve fits associated with body 
line math models. 

Main control routine for computing lift induced drag efficiency "E" 
using a Fourier analysis. 

Solves for the intersection of a line and an ellipse. 

Solves for intersection of two lines. 

Reads all input data except for fuselage math model. Sets up arrays 
and storage areas. Sets up all crude and fine coordinate systems. 

Controls cubic spline fit for interpolation of input spanload distribution. 
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MDOTV 

ORDIN 

PLOTER 

POCRUD 

POFINE 

QWIKDE 

QWIKLO 

REDUCX 

RFAST 

RLORD 

SAREA 

SERIES 

SETNXT 

SIMPSN 

SINCOS 

SLOPBL 

SLOPE 

SLOPY 

SMOTH 

SOLVE B 

SPLINE 

SPLN1 

SPLN1X 

SPRINT 


Performs matrix multiplication of a vector. 

Linear interpolation. 

Controls all graphic output (except input geometry verification). 

Prints results in global crude grid. 

Prints results in wing and body fine grid arrays. Integrates wing 
pressure distributions. Computes body friction drag. 

Main control routine for Quick-Geometry definition and check out. 

Main control routine for interrogation of Quick-Geometry math model. 

Performs interpolation to new grid. 

Fast Fourier transform of real data. 

Bradshaw's L function. 

Computes body surface area given an array of cross-sections. 

Determines Fourier series coefficients. 

Reorders points for non-linear root finder. 

Simpson's rule integration. 

Adjusts input interrogation angles for top and bottom dead center. 

Slope of a tabulated function at an ai'bitrary point. 

Computes boundary conditions for wing surface and axisymmetric bodies. 
Computes wing surface slopes. 

Function for smoothing an array of values. 

Solution of two simultaneous linear algebraic equations. 

Computes a cubic spline through a set of points. 

Computes continuous derivatives interpolation by means of a cubic fit. 
Entry for special cases requiring extrapolation beyond ends of X and Y 
tables. 

Prints output of profile results. 


134 



TANCAL 

Computes characteristic angles for use in the solution (equation 21 of 
Bradshaw and Ferriss) . 

THELIM 

Creates and controls use -theta arrays to establish continuity in the 
cross-sectional model. 

TRID 

Solves tri-diagonal matrix. 

VBRAD 

Computes Bradshaw compressible 2-D turbulent boundary layer 
simulating 3-D boundary layer on infinite yawed wing by Nash-Tseng 
modified chord technique . 

VDOTV 

Computes a vector dot product. 

VINTER 

Performs cubic <5 * fit for separated boundary layer in wing section 
cove regions. 

VNUSUB 

Computes the Nash effective viscosity. 

WINGF 

Finite difference approximations and relaxation solution for wing 
fine grid. 
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KEY VARIABLE DESCRIPTIONS 


This description of key program variables and constants which are located 
several common blocks will be useful in understanding flow logic. 


VARIABLE 

DESCRIPTION 

AK 

The value 1-M . 

A LAM 

Wing taper ratio (A). 

ALPHA 

Angle -of-attack (radians). 

AM AC 

Wing mean aerodynamic chord (MAC). 

AMACH 

Mach number. 

AM2 

2 

The value M . 

AOA 

Angle -of-attack (degrees). 

AR 

Wing aspect ratio (/R). 

BAREA 

Body wetted area. 

BCF 

Body skin friction coefficient. 

BCL 

Wing crude grid lower boundary slopes. 

BCLF 

Wing fine grid lower boundary slopes. 

BCU 

Wing crude grid upper boundary slopes. 

BCUF 

Wing fine grid upper boundary slopes. 

BNOSE 

X-coordinate of body nose. 

BODCD 

Body (integrated) drag coefficient. 

BODCL 

Body (integrated) lift coefficient. 

BODCM 

Body (integrated) moment coefficient. 

BPAREA 

Body projected area. 



VARIABLE 

DESCRIPTION 

BS 

Body plot scaling coefficient. 

BTAIL 

X -coordinate of body tail. 

CA 

Global crude grid stretching coefficient ( ( ). 

CAV 

Wing average chord 

CB 

Global crude grid stretching coefficient (f ). 

CC 

Global crude grid stretching coefficient (7^). 

CD 

Global crude grid stretching coefficient (7yy). 

CDI 

Lift induced drag coefficient. 

CDINT 

Wing section integrated drag. 

CE 

Global crude grid stretching coefficient ( f z )- 

CF 

Global crude grid stretching coefficient ( f ), 

CFINT 

Integrated wing section friction coefficient. 

CIR 

Wing circulation (T). 

CLINT 

Wing section integrated lift. 

CMINT 

Wing section integrated moment. 

CMLOC 

Wing section integrated moment about local quarter chord 

CPL 

Wing lower surface pressure coefficient. 

CPU 

Wing upper surface pressure coefficient. 

CSCUT 

Body x-station for cross-sectional cut. 

DE LSL 

Wing boundary layer slopes for section lower 
surface . 

DELSU 

Wing boundary layer slopes for section upper 
surface . 

DETA 

Global crude grid mesh spacing in 17 direction. 
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VARIABLE 


DESCRIPTION 


DIM 

DRDXC 

DRDXF 

DXB 

DXI 

DXW 

DYB 

DYW 

DZB 

DZETA 

DZW 

E 

ETA 

G 

II 

IBGI 

IBGL 

IL 

ILEF 

IMACH 

IMAXB 

IMAX 

IMAXW 


Configuration length for non-dimensionalizing maximum 
potential updates. 

Axisymmetric body slope distribution in crude grid. 
Axisymmetric body slope distribution in fine grid. 

Fine body grid mesh spacing in X direction. 

Global crude grid mesh spacing in £ direction. 

Fine wing grid mesh spacing in X direction. 

Fine body grid mesh spacing in Y direction. 

Fine wing grid mesh spacing in Y direction. 

Fine body grid mesh spacing in Z direction. 

Global crude grid mesh spacing in f direction. 

Fine wing grid mesh spacing in Z direction. 

Wing spanload efficiency. 

V coordinates for global crude grid (transformed 
space) . 

2 

The value ( 7 + 1) M . 

2 

The value ( 7 -1) M . 

Crude grid I value of body grid inner overlap region (forward). 
Crude grid I value of body grid inner overlap region (aft). 
Crude grid wing leading edge I values. 

Wing fine grid leading edge I value. 

Code for subsonic (0) or supersonic (1) flow at a grid point. 
Maximum number of fine body grid points in X direction. 
Maximum number of crude grid points in X direction. 
Maximum number of fine wing grid points in X direction. 
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VARIABLE 


DESCRIPTION 


INOSE 

INOSEC 

IT 

ITAIL 

ITAILC 

ITEF 

ITER 

JBG 

JMAX 

JMAXB 

JROOT 

JSD 

JSDC 

JTIP 

KBB 

KBOD 

KBC 

KBGL 

KBGU 

KBW 

KLOC 

KLOF 

KM AX 


Fine body grid I value at body nose. 

Crude grid I value at body nose. 

Crude grid wing trailing edge I values. 

Fine body grid I value at body tail. 

Crude grid I value at body tail. 

Wing fine grid trailing edge I value. 

Iteration count. 

Crude grid J value of body grid inner overlap region (side). 
Maximum number of crude grid points in Y direction. 
Maximum number of fine body grid points in Y direction. 
Grid J value at wing root. 

Fine body grid J value at first influence of body boundary 
points. 

Crude grid J value first influenced by body boundary point. 
Grid J value at wing tip. 

Fine body grid K value at wing plane. 

Code for body in crude grid (0) or fine embedded grid (1). 
Crude grid K value at wing plane. 

Crude grid K value of body grid inner overlap region (lower) 
Crude grid K value of body grid inner overlap region (upper) 
Fine wing grid K value at wing plane. 

Crude grid K limiters for body surface (lower). 

Fine grid K limiters for body surface (lower). 

Maximum number of crude grid points in Z direction. 
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VARIABLE 


DESCRIPTION 


KMAXB 

KMAXW 

KODB 

KUPC 

KUPF 

MAXIT 

MAXITF 

MODV 

NCASE 

NINB 

NINW 

NOSEB 

NOSEW 

NPOA 

NPOB 

NSECT 

NTC 

NTF 

NTOTB 


Maximum number of fine body grid points in Z direction. 

Maximum number of fine wing grid points in Z direction. 

Body option code . . . cylinder, axisymmetric , arbitrary 
body (input as BKOD) . 

Crude grid K limiters for body surface (upper). 

Fine grid K limiters for body surface (upper). 

Maximum number of initial crude grid iterations 
(input as AXIT). 

Maximum number of crude/fine grid cycles (input as AXITF) . 

Mode of operation for viscous effects (input as VISMOD). 

Case description . . . wing, body or wing-body case (input 
as CASE). 

Number of ordinates defining axisymmetric body shape 
(input as BNIN) . 

Number of ordinates defining each wing section (input as ANIN). 

Blunt/sharp nose body code (for spline fit). 

Blunt/sharp nose wing code (for spline fit). 

Number of fine grid points between leading and trailing 
edge of each wing section. 

Number of fine body grid points between nose and tail of body. 

Number of defining wing sections (input as ASECT). 

Number of points representing body cross-sections in crude 
grid. 

Number of points representing body cross-sections in fine 
body grid. 

Total number of fine body grid points in single X-Z plane. 
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VARIABLE 


DESCRIPTION 


NTOTAL 

NTOTW 

NWPO 

PBL 

PCI 

PC 2 

PC 3 

PCL 

PF1 

PF2 

PF3 

PFL 

PI 

RADIUS 

RAVC 

RAVF 

RC 

RE 

RF 


Total number of crude grid points in single X-Z plane. 

Total number of fine wing grid points in single X-Z plane. 

Code for print out of crude grid results for diagnostic 
purposes (input as WPO) . 

Fine body grid wing/wake lower surface potentials. 

Global crude grid potential (^) arrays (Note: Only three 
planes are in core at one time.) 

Global crude grid potential (f ) arrays (Note: Only three 
planes are in core at one time.) 

Global crude grid potential (^) arrays (Note: Only three 
planes are in core at one time.) 

Crude grid wing/wake lower surface potentials. 

Fine embedded wing and body potential arrays. (Note: Only 
three planes are in core at one time). 

Fine embedded wing and body potential arrays. (Note: Only 
three planes are in core at one time) . 

Fine embedded wing and body potential arrays. (Note: Only 
three planes are in core at one time) . 

Fine wing grid surface potentials, lower surface 

v 

Radius for body cylinder option. 

Crude grid average body radius for boundary condition calculation. 
Fine grid average body radius for Ixmndary condition calculation. 
Axisymmetric body radius distribution in crude grid. 

Free stream Reynolds number. 

Axisymmetric body radius distribution in fine grid. 
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VARIABLE 


DESCRIPTION 


RIN 

RMAX 

SDD 

SEXP 

SFD 

SGRAD 

SLD 

SMD 

THETC 

THETF 

TITLE 

TSLOC 

TWIST 

W 

WARE A 

WCD 

WCF 

WCL 

WCM 


Input R ordinates defining axisymmetric body. 

Body maximum radius for computational body surface. 
Wing spanwise drag coefficient CC^ . 


C 


AV 


Wing exposed area. 

Wing spanwise friction coefficient CC^ . 


Body side slope at wing -body juncture. 
Wing span load coefficient CC^. 


C 


AV 


Wing spanwise moment coefficient CC 

m 

C AV 

Body crude grid angular cuts. 

Body fine grid angular cuts. 

Case title for identifying graphic and printed output. 

Wing local sweep angle at wing fine grid boundary points. 

Wing twist (incidence) distribution. 

Relaxation factor a> . 

Wing area S, I7 . 

W 

Wing drag coefficient. 

Wing friction drag coefficient. 

Wing lift coefficient. 

Wing moment coefficient. 
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VARIABLE 


DESCRIPTION 


WCORD 

WS 

X 

XBF 

XI 

XI LE 

XINB 

XINW 

XITE 

XLE 

XLET 

XMOM 

XNC 

XNF 

XNOSE 

XOL 

XPL 

XPT 

XSF 

XTAIL 

XTE 

XTET 

XWF 


Wing section local chord length. 

Wing plot scaling coefficient (input as SCALW). 

X coordinate for global crude grid (physical space) . 

X coordinate for body fine grid. 

£ coordinate for global crude grid 
(transformed space) . 

{ coordinate of local wing section leading edge. 

Input X ordinates defining axisymmetric body . 

Wing section defining X-ordinates. 

( coordinate of local wing section trailing edge. 

X coordinate of local wing section leading edge. 

X coordinate of wing tip leading edge. 

Position about which configuration moments are computed. 

Body normal vector X direction at crude grid body points. 

Body normal vector X direction at fine grid body points . 

X coordinate of body nose repositioned in crude grid. 

Non-dimensional distance along body length and wing chord. 

X coordinate of input wing section leading edge. 

X coordinate of input wing section trailing edge. 

X coordinate defining wing sections at each fine grid 
boundary points. 

X coordinate of body tail repositioned in crude grid. 

X coordinate of local wing section trailing edge. 

X coordinate of wing tip trailing edge. 

Fine embedded wing grid X coordinate. 
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variable 


DESCRIPTION 


Y 

YBF 

YINL 

YINU 

YNC 

YNF 

YOB 

YP 

YSF 

YTIP 

Z 

ZBF 

ZETA 

ZNC 

ZNF 

ZWF 

ZWING 


Y coordinate for global crude grid (physical space) . 

Y coordinate for body fine grid* 

Wing section defining lower Y -ordinates. 

Wing section defining upper Y-ordinates. 

Body normal vector Y direction at crude grid body points. 
Body normal vector Y direction at fine grid body points. 
Wing span station (2Y/b). 

Y coordinate of input wing section trailing edge. 

Y coordinate defining wing sections at each fine grid 
boundary points. 

Y coordinate of wing tip. 

Z coordinate for global crude grid (physical space). 

Z coordinates for body fine grid. 

f coordinate for global crude grid (transformed space). 
Body normal vector Z direction at crude grid body points. 

Body normal vector Z direction at fine grid body points. 
Fine embedded wing grid Z coordinate. 

Wing height relative to center of body. 
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RECOMMENDATIONS FOR USAGE 


The thing that becomes apparent when first using a three-dimensional transonic 
relaxation scheme is that a considerable amount of computer time and core is required 
for a solution. Perhaps, an order of magnitude increase over typical resources for a 
subsonic/supersonic panel method will be noted. And, all of this effort will result in 
a solution for a single angle -of-attack, Mach number, and Reynolds number combination. 
It becomes very important to conserve time and resources when using this type of 
methodology. The geometry verification section of the method has been developed for 
this reason. It has been found advantageous to use this code in the following manner. 
When first setting up for analysis of a new configuration, plan on submitting a 
secondary job along with the primary job (which will provide the complete solution). 

The secondary job should be set (AXIT = 0, AXITF = 0) for the geometry verification 
mode or set the solution for a single crude and single fine grid iteration. This will 
provide a complete set up and cycle through the entire code. Since time requirements 
will be low, the secondary job should be returned quickly. If errors are found, the 
primary job can be cancelled to save costs and the error can be quickly corrected. 

Analyses performed on highly swept/highly tapered wing planforms at extreme 
flow conditions indicate that occasionally an erroneous diverging condition at the wing 
tip is possible. This situation can be identified by increasing values of Av«> MAX at the 
last span station on the wing. For severe cases, propagation of the discrepancy in- 
board will be noted. This problem is caused by the large differences in resolution 
between the embedded wing grid system and the surrounding crude grid system. It 
occurs only when strong shock waves exist at the wing tip. A code modification has 
been developed to relieve this problem so solutions at extreme conditions can be 
obtained. In the main program, a parameter "KTIP" has been set to 0. This will 
provide a normal solution process. If wing tip problems are encountered, "KTIP" 
should be set to 1 (KTIP = 1). This will have the effect of reducing wing fine grid 
density in general (80 evenly spaced points along the chord) and severely stepping 
down resolution at the wing tip. This will make the crude and fine grid systems more 
compatible from a numerical standpoint. This option should not be used for the 
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majority of cases for which it will not be required. Note that the included sample cases 
provide examples of geometries and flow conditions for which this option was not 
necessary. 

In most cases, modeling wing -body combinations using the crude grid body 
only option (BKOD < 0) should be quite sufficient. As might be expected, computing 
resources are conserved. There is normally high grid density in the wing -body 
juncture region because of the crude grid stretching with its origin near the juncture. 

If details of the fuselage are important (canopies, blisters, and fairings), the fine 
body grid option is recommended (BKOD > 0). Isolated bodies, because of the reduced 
computing requirements, should always be analyzed using the embedded body grid 
system . 

It is recommended that the wing crude grid system output option (WPO = 1) be 
used. While crude grid wing C p distributions are of little value beyond diagnosing 
errors, the Mach charts will provide a measure of the extent of the supersonic flow 
region into the flow field. 

The wing and wing-body sample cases were all computed by using 100 crude 
iteration cycles followed by 80 crude-fine iteration cycles. Comparisons made with 
experimental data for the sample cases and additional wing-body configurations 
indicated that this level of computational convergence would be satisfactory for 
engineering applications. Any discrepancies between the experiment and analysis 
were attributed to the methods small-disturbance character and the high flow gra- 
dients encountered in certain applications. Waggoner^ 40 \ however, was shown 
that the basic 100/80 iteration cycle count may be insufficient for high aspect ratio 
supercritical wing cases. This may, in part, explain correlation discrepancies 
noted for the NASA supercritical wing-fuselage configuration (page 82) and the 
Douglas wing cylinder configuration (page 92). For this class of geometry, the 
total number of iterations may have to be increased by a factor of 3 or 4. The 
user should be aware of this lack of convergence if applications to high aspect 
ratio wings or supercritical type sections are of interest. 
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